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Real Gas Effects in Flow Over Blunt Bodies 
at Hypersonic Speeds' 


H. T. NAGAMATSU,* R. E. GEIGER,** ann R. E. SHEER, JR.*** 
General Electric Company 


Summary 


A hypersonic shock tunnel has been developed to investigate 
the aerodynamie characteristics of flow over bodies at conditions 
comparable to those encountered by ballistic missiles and satel- 
lites re-entering the atmosphere. Some results for a shock veloc- 
ity of over 50,000 ft. per sec. in the shock tube portion of the 
facility are presented. Static pressure investigations were made 
in the nozzle to determine the flow condition and the expansion 
process. 

The results of the investigation of representative blunt bodies 
at hypersonic Mach numbers and nozzle stagnation temperatures 
up to approximately 6,000°K. are presented. These include 
body pressure distributions, shock-wave shapes, detachment dis- 
tances, and photographs of the luminous gas region in the shock 
layer. It is seen that the shock detachment distance is smaller 
at higher stagnation temperatures due to the real gas effects. 
For the hemisphere the pressure distribution was less than that 
predicted by the modified Newtonian theory for all stagnation 
temperatures. For a 50° cone-hemisphere the pressure distribu- 
tion and the shock detachment distance were appreciably 
affected by the real gas effects. 

The observed shock-wave shape and the approximate boundary 
layer on a flat plate are compared with the analytical prediction. 
Some preliminary results for the detached shock wave produced 
by a blunt two-dimensional body in a low density flow at a Mach 
number of 19.6 are presented. 


(1) Introduction 


: problem of obtaining aerodynamic information 
at hypersonic flight speeds has currently received 
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1959. Revised and received July 1, 1959. 
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considerable interest in connection with the develop- 
ment of intercontinental ballistic missiles, satellites, 
and space vehicles. For flight Mach numbers greater 
than about 10, air can no longer be considered as 
being a simple mixture of diatomic oxygen and nitro- 
gen. At a Mach number of 20 the air temperature 
behind a normal shock wave for a blunt body can be 
as high as 6,500°K. during re-entry. At these high 
temperatures the effects of dissociation, ionization, 
and chemical reaction become very great. These high- 
temperature phenomena are referred to as ‘‘real gas 
effects."” Knowledge regarding the real gas effects 
upon the flow over a body is important for the design 
of efficient hypersonic vehicles. 

One test facility which can simulate hypersonic flight 
Mach numbers and temperatures is the hypersonic 
shock tunnel, one of which was completed at the General 
Electric Research Laboratory in 1957. Initial opera- 
tion of this facility was devoted to an extensive series 
of tests designed primarily to explore the facility and 
instrumentation performance and to establish operating 
procedures. More recently, investigations have been 
conducted over a range of temperatures and flow Mach 
numbers with various models in the test section. 

The purpose of this paper is to present some of the 
recent results obtained in the hypersonic shock tunnel. 
Some investigations of the state of the gas during the 
expansion in the nozzle will be discussed. The real gas 
effects at high temperatures upon the detached shock- 
wave shape and pressure distribution for a hemi- 
sphere and a 50° cone-hemisphere will be presented. 
Finally, some results for nearly perfect gas flows over a 
flat plate and various blunt bodies at hypersonic Mach 
numbers will be discussed. 
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Fic. 1. Idealized picture of shock tunnel phenomena. 


(2) Description of Facility, Instrumentation, 
Models, and Test Conditions 


(a) Hypersonic Shock Tunnel 


A shock tube!~ consists of a high-pressure or driver 
section and a low-pressure or driven section separated 
by a diaphragm. If a nozzle is placed at the end of the 
shock tube, the compressed and heated gas behind the 
incident wave may be expanded to high velocities in a 
nozzle to produce hypersonic flow in the test section. 
This arrangement is referred to as a shock tunnel. If 
the flow behind the incident shock wave is passed 
directly into a diverging section, this arrangement is 
known as a nonreflected or straight-through nozzle.” * 
If the shock tube is closed except for a small entrance to 
a convergent-divergent nozzle, the incident shock wave 
reflects from the end of the tube, which further com- 
presses the gas and brings the flow almost to rest. 
This doubly heated and compressed gas then expands in 
the nozzle to produce hypersonic flow at high stagna- 
tion temperatures. Such an arrangement is called a 
reflected nozzle.* A schematic drawing of the re- 
flected shock tunnel processes is presented in Fig. 1. 

The Research Laboratory shock tunnel (Fig. 2) 
consists of a 20-ft. driver, a 103-ft. driven tube, a 4-ft. 
reflected type conical nozzle, and a 5-ft. diameter verti- 
cal dump tank. A more detailed description of this 
facility is presented in reference 10. All shock tunnel 
experiments to date have used combustion driving. 
Stoichiometric mixtures of hydrogen and oxygen diluted 
by helium!'! have been used with ignition accomplished 
by eighteen spark plugs distributed along the driver 
walls in a spiral array. 


(b) Instrumentation 


Quartz and barium-titanate piezoelectric transducers 
are used to measure pressures on the shock tube walls, 
in the static pressure probe, and on the models in the 
test section. 

These gages are calibrated in the models directly be- 
fore and after each series of tests by mounting the 
models at the end of the shock tube and calibrating with 
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Fic. 2. General Electric Research hypersonic shock 
tunnel. 


the pressure jump across weak reflected shocks produced 
by an incident shock Mach number of about 1.4. 

Thin platinum film heat gages, of the type described 
in reference 12-15, are used to qualitatively observe 
wall temperature histories and as time of arrival indi- 
cators in shock velocity measurements. 


(c) Models 

The hemisphere pressure distribution data were 
obtained with a 5-in. diameter brass model with pres- 
sure orifices located at spherical ray angles of 0°, 30°, 
45° (2), 60°, and 80°. 

Two brass 5-in. diameter models were used to meas- 
ure the pressure on the 50° cone-sphere configuration. 
One model with gages located at S/R = 0, 0.145, 
0.522, and 1.63, and the other with S/R = 0, 0.45, 
0.85, and 1.40, where S is the distance along the model 
surface from the stagnation point and R is the spherical 
nose radius. 

A static pressure probe with a pointed nose section 
and piezoelectric pickup was used to measure static 
pressure in the nozzle. 


d) Test Conditions 


The conditions in the test section are essentially 
determined by the nozzle area ratio and conditions 
behind the reflected shock at the downstream end of the 
shock tube. The strength of the shock wave is con- 
trolled by the choice of the driver loading and the 
initial pressure in the driven part of the shock tube. 
The initial temperature in the driven tube is always 
room temperature. 

The conditions behind the reflected shock wave are 
the stagnation conditions for the nozzle expansion 
process. Pressures behind the reflected shock wave 
are measured directly with the standard Kistler quartz 
pickup, whereas the reflected temperatures, 7;, are 


calculated values, obtained by numerically solving the 
normal shock equation for both incident and reflected 
shocks, assuming the gas to be in thermodynamic 
equilibrium. 


The thermodynamic data of references 
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16-19 are used in the calculations. The shock-wave 
velocity along the driven tube is measured at twelve 
stations by using a modified Tektronix oscilloscope, 
and at the end of the tube the shock velocity is meas- 
ured by a Berkeley counter. This velocity is used to 
calculate the reflected equilibrium stagnation tempera- 
ture for the nozzle. With the combustion driver 
technique, the shock-wave attenuation in the down- 
stream portion of the driven tube is much less than 
that occurring near the beginning of the tube. Aill 
such temperatures referred to in this paper were de- 
termined in this manner. 

The majority of investigations to date have been con- 
ducted with a 1-in. diameter throat section to give a 
geometric area ratio of 576. The corresponding flow 
Mach number at the nozzle exit according to one-dimen- 
sional ideal inviscid gas theory is 10.15 for y = 1.4. 
This configuration is referred to as the nominal Mach 
10 nozzle. Some exploratory investigations have been 
conducted with nominal Mach 16, 20, and 26 nozzles. 

For relatively low temperatures, where air still 
behaves as an ideal gas, test section Mach numbers 
are determined by measuring the model impact pres- 
sures. At elevated temperatures the flow Mach num- 
ber is estimated by assuming thermodynamic equilil» 
rium and using the results of a series of one-dimen- 
sional isentropic nozzle calculations for 1,000 psia 
stagnation pressures covering a range of temperatures 
up to 7,000°K. These calculations used the air data 
of references 16-20). At high reflected stagnation pres- 
sures the flow in the nozzle was observed to expand 
under nearly thermodynamic equilibrium conditions 
as discussed in the following section. 


(3) Flow in Nozzle 
(a) Flow Uniformity in the Test Section 


By operating the shock tunnel such that the re- 
flected temperature at the entrance to the nozzle is 
1,000 to 1,500°K. and the reflected pressure is 500 to 
2,000 psia, the real gas effects are very small. The 
flow condition in the test section was investigated by 
measuring the impact pressure distribution, the shock- 
wave shape for a two-dimensional 40° half-angle wedge 
and a 4-in. diameter 60° cone, and the static pressure 
in the test section. 

The impact pressure distribution over a 5-in. span 
indicated a fairly uniform flow." A centerline impact 
Mach number of 10.0 was determined from the meas- 
ured impact pressure and the reflected pressure for a 
reflected stagnation temperature of 1,100°K. By 
measuring the centerline static pressure in the test 
section, the Mach number determined from the ratio 
of static to reflected pressure is 10.1. Thus, the Mach 
number determined by the two methods agrees very 
closely for the condition of nearly perfect gas flow 
conditions. 

The shock waves for the two-dimensional wedge and 
the 60° cone were nearly straight up to the point of 
interaction with the expansion wave emanating from the 


shoulder. This indicates that there were no large dis- 
turbances in the test section. 

The impact pressures indicated that the test section 
Mach number is nearly constant for 1 to 4 millisec., 
depending upon the operating conditions. It requires 
usually 0.5 to 0.7 millisec. to establish flow in the test 
section and, consequently, there is adequate time for 
obtaining good measurements. 


(b) Nozzle Static Pressure and Wedge Shock Angle 
at High Temperatures 

For the model results to be meaningful, the state of 
the air in the test section must be known. At high 
reflected temperatures, the air dissociates, ionizes, and 
chemical reactions take place between the constituents. 
As the heated air is expanded in the nozzle, the gas can 
expand in thermodynamic equilibrium, or, if the den- 
sities are low enough, the expansion may be too rapid 
to permit the attainment of an equilibrium gas state. 

In order to ascertain the magnitude of the real gas 
effects on the nozzle expansion, calculations were made 
for two extreme cases—equilibrium nozzle and model 
flow, and completely “‘frozen’’* nozzle flow with equi- 
librium over the model. These calculations assumed 
one-dimensional flow with a sonic throat. The 
equilibrium properties of air were obtained from ref- 
erences 17-20). Some of the computational results are 
presented in Fig. 3. They show that for a given nozzle 
area ratio, the model impact pressure is not too sensi- 
tive to whether the nozzle expansion process is an 
equilibrium or frozen one. However, nozzle static 
pressure and two-dimensional wedge shock angles are 
calculated to be quite different for the two extremes of 
equilibrium and completely frozen nozzle expansions. 
Accordingly, a number of experiments were performed 
in which measurements were made of the centerline 
static pressure at a location halfway into the nominal 


* As used in this paper, a completely frozen expansion is one in 
which only the rotational and translational energies of the gas 
adjust during the expansion while the energy content cf all other 
degrees of freedom remain at their stagnation values. 
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Fic. 3. Isentropic pressure ratio in nozzle vs. geometric nozzle 
area ratio. 
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Fic. 4. 
' temperature at area ratio of 144. 


Mach 10 nozzle, and of the shock wave angle for a 22° 
half-angle wedge located 1 in. downstream of the noz- 
zle exit. 


The static pressure results are presented in Fig. 4. 
Also shown are the calculated isentropic one-dimen- 
sional flow values for the cases of equilibrium and com- 
pletely frozen expansions from the nozzle stagnation 
conditions—.e., from equilibrium conditions behind the 
reflected shock wave. These calculated curves do not 
contain an estimate of the effect of changing boundary- 
layer thickness. The veritical bars through the data 
points represent the variation with time observed during 
the individual experiments. 


The nozzle expansion calculations indicated that the 
ratio of static pressure to stagnation pressure is con- 
siderably higher if the gas is in thermodynamic 
equilibrium rather than if it is in frozen flow. The 
static pressure results, Fig. 4, indicate that at a given 
reflected temperature higher ratios of static to stagna- 
tion pressure are observed as the reflected pressure is 
increased, up to a point, but further increases have no 
effect. At lower driven pressure, which corresponds to 
a lower reflected pressure, the static to reflected pressure 
ratio approaches the frozen flow curve for a given re- 
flected temperature, Fig. 4. As the reflected tempera- 
ture is increased, the reflected pressure for which 
nearly equilibrium expansion takes place in the nozzle 
becomes higher. 


With 400 psia charging pressure in the driver, the 
nozzle flow appears to be nearly equilibrium for re- 
flected temperatures up to 5,200°K. At higher 
temperatures the measured pressures in the nozzle 
approached the frozen expansion values. Further 
experiments will be conducted at higher reflected 
temperatures and higher driver pressures to obtain 
more knowledge of the state of gas during the expan- 
sion process. 


The wedge shock angles were measured at about 
5,000°K. reflected temperatures with driver loadings 
of 400 and 240 psia. The observed shock wave angles 
agreed with the calculated curve for equilibrium during 
the nozzle expansion and also after the shock wave. 
The measured shock angles for the two different 
driver pressures show no difference within the accuracy 
of the measurement. An approximate correction for 
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the wedge boundary layer was determined from the 
schlieren photographs. 

The shock angles observed on the 22° wedge rules out 
the possibility of a frozen nozzle expansion with equi- 
librium behind the wedge shock. The nozzle static 
pressure results, Fig. 4, indicate that the expansion in 
the nozzle takes place very nearly at equilibrium for the 
conditions of these wedge experiments. 

Thus, the preliminary results of the static pressure 
and wedge shock wave angle measurements appear to 
indicate that the expansion of the air in the nozzle is 
essentially in equilibrium for the higher operating pres- 
sure at the various reflected temperatures up to about 
5,200°K. In general, as the temperature is increased, 
higher operating pressures are required to attain an 
equilibrium or nearly equilibrium expansion. For low 
reflected pressures at a high temperature, the meas- 
ured static pressures approach the calculated values 
for a frozen expansion. ‘Therefore, the pressure distri- 
bution and schlieren photographs for models in the 
test section were obtained at high enough operating 
pressures to ensure nearly equilibrium flow in the 
nozzle. 

In Table 1 the flow Mach number, MM, ambient 
temperature, 7’, static pressure, p, density ratio across 
a normal shock wave, p2/p, and stagnation temperature 
after a normal shock wave, 7:a,, are presented as a 
function of the reflected equilibrium stagnation 
temperature, 7, at the entrance to the nozzle. 


(4) Real Gas Effects on Hemisphere Shock Shape 
and Surface Pressure 


A series of tests were conducted with a hemisphere 
model in the test section of a nominal Mach 10) nozzle. 
The calculated reflected stagnation temperature at the 
entrance to the nozzle was varied between 1,000 and 
6,000°K. Schlieren photographs and photographs 
of the luminosity of the gas itself as well as the measure- 
ments of the surface pressures were obtained. 

Figs. 5(a) and 5(b) are schlieren photographs of the 
detached shock waves for reflected stagnation tempera- 
tures of 1,060 and 4,290°K. In Fig. 5(b) the lumino- 
sity of the hot gas back of the detached shock wave can 
be observed. At this temperature there is a short 
dark region between the shock and the beginning of the 
luminosity, suggesting relaxation of the excited species 
producing the visible radiation. 

In Fig. 6 the two shock waves are compared with 
each other and with that obtained at a Mach number 
of 5.8 under ideal gas flow conditions.*!_ Two interest- 


ing features are indicated in this figure. First, the 
TABLE | 
7; M Ppsia po/p 
1,100 10 59 0.0134 5.70 1,100 
2,000 9.4 115 0.0214 6.20 1,960 
3,000 8.4 244 0.0170 7.50 2,850 
4,000 7.8 450 0.0206 10.2 3,370 
4,290 7.6 480 0.0227 10.5 3,600 
5,300 6.1 1310 0.0114 11.3 4,580 


| 
| 


Fic. 5( 


t 
Fic. 
f 
| 
' 
Be 


the 


Ps Out 
equi- 
static 
on in 
or the 


bssure 
par to 
zle is 
pres- 
bout 
based, 
in an 
r low 
neas- 
alues 
listri- 
1 the 
ating 
the 


bient 
LCTOSS 
ature 

as a 
ation 


hape 


phere 
ozzle. 
t the 
) and 
raphs 
sure- 


f the 
pera- 
nino- 
e can 
short 
f the 


with 
m ber 
‘rest- 
the 


REAL GAS EFFECTS 245 


HEMISPHERE SHOCK SHAPE 


1060°K, M=10 
“154290 °K, Meg 7.5(CALC) 
CALC. SONIC 
POINTS / 
/ 
/ 


——— PRESENT RESULTS 
21, M=5.8 


| 
| 
| 


Fic. 6. Detached shock wave shapes for hemisphere at different 
temperatures. 
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Fic. 5(a). Hemisphere schlieren picture at 1,060°K. reflected 
stagnation temperature. 
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Fic. 5(b). Hemisphere schlieren picture at 4,290°K. reflected Fic. 8. Luminosity photograph of hemisphere at 5,800°K. 
stagnation temperature. reflected stagnation temperature. 
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FIG. 


detached shock shapes are very similar despite the 
difference in Mach number and stagnation temperature. 
Secondly, a comparison of the detachment distance 
for the 1,060 and 4,290°K. shock waves, where the 
corresponding equilibrium flow Mach numbers are 10 
and 7.6, respectively, indicates that the higher tempera- 
ture effect, which causes a decrease in detachment dis- 
tance, dominates the lower Mach number effect, which 
causes an increase in detachment distance. This is 
in agreement with the theoretical results which indicate 
that the density ratio across the shock wave is the 
primary parameter in determining the blunt body shock 
detachment distance. The density ratio across a nor- 
mal shock wave for an ideal gas varies very little with 
Mach number in the range above 8, whereas for air at a 
given Mach number it increases appreciably with in- 
creasing stagnation temperature (Fig. 7). 

A photograph showing the luminosity of the gas was 
taken at a reflected stagnation temperature of 5,S00°K. 
(Fig. 8). The luminous region seems to extend to the 
shock wave, and the dark region after the shock wave 
observed at 4,290°K. does not seem to be present. 
The indication of decreased relaxation time at higher 
temperatures after the normal shock wave has been 
observed after the shock wave in a constant area 
tube.”* °° As the gas expands around the body from 
the stagnation point and the detached shock wave be- 
comes weaker, the gas cools and loses its luminosity. 

The shock detachment distances obtained from the 
schlieren photographs taken at 1,060°K. and 4,290°K. 
are compared with some theoretical values in Table 
2. Thermodynamic equilibrium was assumed in cal- 
culating the shock density ratios used in the theoretical 
equations for the 4,290°K. case. 

In Fig. 9 the detachment distances based upon 
reference 27 for the cases of equilibrium in the nozzle 
and after the normal shock wave and for completely 
frozen flow in the nozzle with equilibrium after the 
detached shock wave are compared with experimental 
values measured from schlieren and _ self-illuminated 
photographs over a range of temperatures. It is ap- 
parent that as the reflected stagnation temperature 
increases, the detachment distance decreases because of 
the increase in the shock density ratio. At the begin- 
ning of visible luminosity, about 4,000°K. for these 
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TABLE 2 
T; = T: 
Theory Reference 1,060°K. 4290°K. 

Heybey 24 0.123 ~ 
Li-Geiger 25 0.113 0.073 
Van Dyke 26 0.135 
Serbin 27 0.148 0.079 
Experiment 0.132 0.110 


particular test conditions, there is a slight delay be- 
hind the shock wave before the gas shows this visible 
luminosity, and the front of the luminous region is not 
well defined. At temperatures above about 5,000°K. 
where the shock itself cannot be seen, the fronts of the 
luminous gas region become sharp and the values de- 
termined from the direct luminosity photographs fall 
more in line with those obtained from the schlieren 
photographs at lower temperatures. All the observed 
detachment distances in the present temperature range 
fall between the values predicted by the theory” for 
equilibrium and frozen flow in the nozzle. The ob- 
served variation of the detachment distance with in- 
creased stagnation temperature is less than that pre- 
dicted by the simplified detachment theories which 
satisfies only the conditions on the axis. To obtain 
better correlation with the experimental data, it will 
be necessary to solve the complete flow field*® °° be- 
tween the shock wave and the body including the 
correct thermodynamic properties of the air. 

Pressures on the hemisphere were measured at 30°, 
45°, 60°, and 8)° from the stagnation point for various 
reflected stagnation temperatures. In Figs. 10(2) 
and 10(b) the measured pressure coefficients are plotted 
as functions of the angle from the stagnation point at 
different temperatures. It is evident from these 
figures that the pressure coefficient for all the locations 
except the 80° point are less than the modified New- 
tonian* at all temperatures and that there is no sig- 
nificant temperature trend. Even at a nozzle stagna- 
tion temperature of 5,300°K. where the gas behavior is 
highly nonideal, as evidenced by a calculated equilib- 
rium normal shock density ratio of 11.3, the pressure 
coefficient distribution is very close to that for the 
nearly ideal gas case of 1,400°K. The observed experi- 
mental results showing the pressure coefficient over 
most of the hemisphere to be less than the modified 
Newtonian values are consistent with recent analytical 
work.” 28, 29 

The pressure behind the observed shock wave at 
1,060°K., Fig. 5(a), has been calculated by assuming 
y = 1.4 and is compared with experimental curves of 
surface pressure distributions in Fig. 11. The radial 
pressure gradient in the shock layer is seen to be very 
small in the vicinity of the calculated sonic point on the 
detached shock wave. It is also evident in this figure 
that the sonic point on the shock wave is much closer 
to the stagnation point than the body sonic point 
determined from the pressure distribution. 


* Private communication from Dr. Eva Winkler from the 
Naval Ordnance Laboratory, White Oak, Md., has indicated 
similar results in the hypersonic wind tunnel for the hemisphere. 
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50° cone-hemisphere schlieren picture at 1,300°K. 
reflected stagnation temperature. 


Fic. 12(a). 


Fic. 12(b). 50° cone-hemisphere schlieren picture at 4,000°K. 
reflected stagnation temperature. 
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Fic. 13. Variation in shock detachment distance for 50° cone- 
hemisphere with temperature. 


It has been observed that the detached shock-wave 
shape for the hemisphere does not change significantly 
within the present range of stagnation temperatures 
but moves closer to the body as the stagnation tempera- 
ture is increased. Thus, the local shock-wave angle 
determined from Fig. 5(a) for the nearly ideal flow 
condition, 7; = 1,060°K., should be a good approxima- 
tion for the entire temperature range of the tests. 
For given free-stream conditions, the pressure behind 
an oblique shock wave varies primarily as the sin® 2, 
where @ is the shock-wave angle. As indicated in Fig. 
11 the pressure gradient across the shock layer is not 
great; hence, it might be expected that at higher tem- 
peratures the hemisphere pressure distribution would 
not vary significantly from the ideal gas case. This is 
indeed what has been observed in these shock tun- 
nel tests at reflected stagnation temperatures up to 
5,300° K. 


(5) Real Gas Effects on 50° Cone-Hemisphere 
Shock Shape and Surface Pressure 


In addition to the hemisphere, a 50° cone-hemisphere 
with a ratio of nose to base radii of 0.623 was investi- 
gated over a range of reflected stagnation temperatures 
in the nominal Mach 10 nozzle. 

Schlieren photographs of the shock wave produced 
by this body are shown in Figs. 12(a) and 12(b) for 
reflected stagnation temperatures of 1,300 and 4,000°K. 
In the latter figure the luminous gas region can again 
be seen. At 1,300°K. the detached shock wave has an 
inflection point over the conical portion of the model; 
similar results have been observed in reference 9 at 
slightly higher Mach numbers. The shock shapes at 
these two temperatures are quite similar, but at the 
higher temperatures the inflection point is not as pro- 
nounced and the shock wave is closer to the body. 

Measurements of the shock detachment distance and 
the luminous gas fronts made from the schlieren and 
luminosity photographs obtained over a range of flow 
reflected stagnation temperatures are presented in Fig. 
13. It is seen that as the temperature was increased 
the shock moved closer to the body due to the increase 
in density ratio across the shock wave. The experi- 
mental values fall between the theoretical values calcu- 
lated by using the results of reference 27 for the cases 
of equilibrium and frozen nozzle expansions. 


The effect of temperature on the surface pressure dis- 
tribution, which compares the measured values of the 
pressure coefficient ratio, C,/C,,,,,, with the modified 
Newtonian values is presented in Figs. 14(a) and 14(b) 
for reflected stagnation temperatures of 1,400, 3,000, 
4,300, and 5,300°K. The solid curve in these figures is 
the modified Newtonian pressure distribution and the 
dashed curve is the experimental result obtained for the 
hemisphere, Fig. 10. At 1,400°K. the pressure coeffi- 
cients are close to the modified Newtonian theory over 
most of the hemispherical part, Fig. 14(a), and are 
above the theory over most of the conical surface. 
Similar results were obtained for a 40° cone-hemisphere 
at a Mach number of 5.8 in a hypersonic wind tunnel.*! 
As the temperature is increased to 3,000°K., the pres- 
sure over the hemispherical part shifts below the modi- 
fied Newtonian value and approaches the results ob- 
tained for a hemisphere over a range of temperatures. 
The pressure over the conical part approaches the New- 
tonian result. This change in the pressure distribution 
at elevated temperatures is due to the real gas effects 
which cause the shock wave to move closer to the body 
and the effective ratio of specific heats to decrease. 
At 4,300°K. the result is very similar to the 3,000°K. 
case. 

For 5,300°K. the pressure on the hemispherical part 
is again below the modified Newtonian, Fig. 14(b), 
but agrees with the data for the hemisphere at the 
same temperature. On the conical surface the pres- 
sure is slightly below the modified Newtonian values. 
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Fic. 14(a). 50° cone-hemisphere pressure distribution for 
reflected stagnation temperatures of 1,400°K. and 3,000°K. 
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Fic. 14(b). 50° cone-hemisphere pressure distribution for re- 
flected stagnation temperatures of 4,300°K. and 5,300°K. 
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TABLE 3 
+ 1.15 1.20 1.30 1.40 
0.574 0.564 0.546 0. F28 
72.4° 69.6° 67.2° 


The detached shock wave for this temperature, as 
indicated by the front of the luminous gas region, 
is closest to the body and the value of 6/R was approxi- 
mately the same for the hemisphere at this temperature, 
Fig. 9. 


(6) Real Gas Effects on Sonic Location 


One of the features of interest in the study of blunt 
body flow fields at hypersonic Mach numbers is the 
location of the sonic point, both on the detached shock 
wave and on the body. Fora gas with constant specific 
heat ratio the shock angle at the sonic point is easily 
calculated and for the low-temperature hemisphere 
tests at a Mach number of 10 it is 67.3°. This occurs 
on the shock wave at a spherical angle of 26.8° from 
the stagnation point. The body sonic location is at 
about 41.7° as indicated by the experimental body pres- 
sure distribution; thus, it is considerably further out 
from the stagnation streamline than the shock sonic 
point. The portion of the detached shock wave ad- 
jacent to the spherical nose of the 50° cone-hemisphere 
is quite similar in shape to that produced by the hemi- 
sphere and its sonic point occurred at about 29.2° for the 
low-temperature tests. The body sonic point for this 
configuration is at the base of the conical surface. 

At higher temperatures where the gas behavior is no 
longer ideal, the shock angle at the sonic point must 
be determined by a separate numerical solution of the 
shock equations for each nozzle flow stagnation con- 
dition. One such calculation was performed for a 
reflected stagnation temperature of 5,000°K. and corre- 
sponding pressure of 1,000 psia. The shock angle at the 
sonic point was calculated to be 72.3° by assuming 
thermodynamic equilibrium to exist in the nozzle 
expansion and behind the detached shock. Since for 
both the hemisphere and 50° cone-hemisphere the de- 
tached shock-wave shapes changed very slightly up to 
about 4,000°K., and since the nozzle flow investigations 
indicate that thermodynamic equilibrium is nearly 
achieved at high reflected pressures, it appears that the 
shock sonic location tends to move toward the axis for 
these two bodies as the flow stagnation temperature in- 
creases. 

The critical or sonic pressure ratio on the surface of 
a blunt body in a high-temperature flow may be ap- 
proximated by using the ideal gas equations with a 
suitable effective value of specific heat ratio. In Table 
3 the values of sonic pressure ratio for different 
values of y are presented together with the shock-wave 
angle at the sonic point. It is seen that the critical 
pressure ratio increases with decreasing ratio of specific 
heats. For air the effective specific heat ratio de- 
creases as the flow stagnation temperature increases. 
Thus, since the measured ratio of pressure coefficients, 


Fic, 15. Flow over a flat plate with sharp leading edge at flow 
Mach number of 10 with 7; = 1,500°K. 


C,/C>,,.. for a hemispkere changes very little with 
temperature, the hemisphere body sonic point tends to 
move toward the stagnation point just as does the shock 
sonic point. In the case of the 5° cone-hemisphere 
body the experimental body pressure ratios indicate 
that the local Mach number on the conical surface 
remains subsonic within the temperature range investi- 
gated. 


(7) Flow Over a Flat Plate and Ultra-High-Speed 
Flows in a Shock Tube and Shock Tunnel 


(a) Flow Over Flat Plate 


A flat plate with a 5-in. span was tested in the Mach 
10 nozzle at a stagnation temperature of 1,500°K. 
and a Reynolds number per inch based upon the free- 
stream conditions in the test section of 60,000. The 
interaction of the hypersonic boundary layer with the 
leading-edge shock wave can be seen in Fig. 15, which 
clearly shows the sharp density, and, hence, temperature 
gradient, existing in the outer part of the thick laminar 
boundary layer. The boundary-layer thickness as 
approximately indicated by the sharp density gradient 
and the shock shape vary as X°-"4 for 0.17 < X¥ < 1.0 
and 0.17 < X < 0.70 in., respectively, and are not 
discernible closer to the leading edge. For this particu- 
lar test condition the surface temperature was close to 
room temperature; thus, it was relatively cool compared 
to the recovery temperature. Even at these conditions 
the laminar boundary layer is very thick. The 0.74 
power boundary-layer growth in the strong interaction 
region is in good agreement with the original theoretical 
result of reference 30. The case of a noninsulated flat 
plate for hypersonic flow has been developed in refer- 
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Fic. 16. Schlieren picture of right circular cylinder at flow Mach 
number of 19.6. 


erence 31, where a boundary-layer growth proportional 
to X*/4 was assumed. At higher flow Mach numbers 
the theory presented in reference 30 indicates that for a 
given Reynolds number the interaction of the shock 
wave and the boundary layer extends to a greater dis- 
tance from the leading edge of the flat plate. 


(b) Ultra-High-Speed Shock Wave in the Shock Tube 


In the driven section of the shock tube, a maximum 
shock velocity to date of 55,000 ft. per sec. has been 
produced in air with the combustion driver technique. 
The induced velocity after the shock wave, assuming 
equilibrium, was calculated to be 52,000 ft. per sec. 
and the corresponding temperature 16,400°K. For this 
particular case the shock wave required approximately 
10 ft. to attain maximum velocity’ * and for 60 ft. the 
shock moved with very little attenuation. 


A discussion of the method of measuring and pro- 
ducing these high-speed shock waves is presented in 
reference 10. The velocity of these shock waves are 
about twice the values predicted by the simplified 
shock tube theory. The effects of the diaphragm 
opening time, mixing at the contact zone, and nearly 
constant pressure combustion are not included in the 
above theory. All of these effects tend to produce 
stronger shock waves, but the last effect, produced by 
letting the diaphragm break at approximately one-half 
of the maximum pressure at the end of burning, gives 
the largest increase in the shock velocity. The shock 
velocity along the 100-ft. tube was measured at twelve 
stations. With our combustion driver technique, it 
has been possible to obtain reproducible strong shock 
waves. 
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It is anticipated that future operation at higher driver 
pressures will permit the attainment of even greater 
shock velocities. Thus, the constant area section of the 
shock tube can be used to obtain stagnation heat-trans- 
fer information corresponding to flight Mach numbers 
of about 50 at high altitude conditions. Moreover, the 
shock tunnel could be filled with gases corresponding to 
the planetary atmospheres of Venus and Mars to ob- 
tain heat-transfer and aerodynamic information for 
those entry problems. 


(c) Preliminary Results at Flow Mach Numbers Above 10 


Some preliminary tests have been made in the higher 
Mach number nozzles with reflected equilibrium 
temperatures of 1,100 and 2,200°K. and pressures from 
1000—2,000 psia where air may still be considered al- 
most an ideal gas. These tests were conducted to as- 
certain whether full flow could be established in the 
nominal Mach 16, 20, and 26 nozzles. The stagnation 
point pressures measured on the various blunt models 
used in these tests indicated no difficulty in establishing 
the flow, provided the dump tank pressure was below 
a certain limiting value. 


For the nominal Mach 20 nozzle the throat diam- 
eter is only 0.190 in. with a nozzle exit diameter of 
24 in. Blunt body impact pressures measured in this 
nozzle indicated the establishment of flow in about 500 
microsec. The flow Mach number based on the impact 
pressure was approximately 19.6. The calculated 
static pressure in the test section assuming isentropic ex- 
pansion was 14 microns of Hg. At this pressure the flow 
in the test section corresponds to conditions encoun- 
tered by vehicles at high altitudes. It was possible to 
observe a detached shock wave for a right circular cylin- 
der of 2-in. diameter, Fig. 16. Thus, at high Mach 
numbers it may be possible to have a rarefied gas flow 
condition ahead of the detached shock wave and a con- 
tinuum condition after the shock.** At lower Mach 
numbers for rarefied flows the shock wave is not sharp 
but covers some finite distance as discussed in refer- 
ence 34. 


(8) Conclusions 


Preliminary results have been obtained in a shock 
tunnel in which it is possible to produce high flow Mach 
numbers with the corresponding high temperatures 
encountered in flight. 

There is no difficulty in establishing hypersonic flow 
in the nozzle as long as the nozzle and the dump tank 
are evacuated to sufficiently low pressures. The dura- 
tion of the steady flow, as determined by the model 
stagnation pressure, ranges approximately from 1 to 
4 millisec., depending upon the stagnation temperature 
and the Mach number. 

The expansion of the heated air in the nozzle appears 
to follow thermal equilibrium for high reflected pres- 
sures, but at lower pressures the expanding gas tends to 
approach the frozen flow condition. 

For a hemisphere in the nominal Mach 10 nozzle the 
real gas effects are to move the shock wave toward the 
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body without any appreciable change in the shape. 
The surface pressure is less than modified Newtonian 


value up to a spherical angle of about 80°. The pres- 
sure distribution remains nearly the same over a range 
of stagnation temperatures up to 5,300°K. 

With the 50° cone-hemisphere much greater real gas 
effects are indicated upon the pressure distribution than 
shown with the hemisphere. The detached shock wave 
moves closer to the body with increasing stagnation 
temperature. For 1,400°K. reflected equilibrium tem- 
perature, the pressure on the conical surface is higher 
than predicted by modified Newtonian theory and the 
pressure on the hemisphere is close to Newtonian. At 
higher stagnation temperatures the pressure on the 
spherical part approaches the hemisphere value and the 
pressure on the cone approaches the modified Newtonian 
value. Thus, for this type of blunt body the high- 
temperature effects upon the flow field appear to be 
appreciable. 

Observation of flow over a flat plate for a nearly 
ideal gas condition showed that a strong shock-wave 
boundary-layer interaction exists at a Mach number of 
10 and that the laminar boundary layer is still very 
thick even with surface temperatures considerably 
less than the recovery temperature. 

An ultra-high-speed shock velocity of 55,000 ft. per 
sec. has been produced in the driven tube with a calcu- 
lated air velocity after the shock wave of 52,000 ft. per 
sec. Thus, stagnation-point heat-transfer rates and 
temperatures encountered by a space vehicle at re- 
entry at high altitudes can be produced in the shock- 
tube portion of the shock tunnel. 

A flow Mach number of almost 20 has been achieved 
in the test section of the nozzle. Even higher flow Mach 
numbers may be possible at high-altitude rarefied 
conditions. 
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Flow About an Unsteadily Rotating Disc 


E. M. SPARROW* anp J. L. GREGG** 
Lewts Research Center, NASA 


Summary 


An analysis is made of the unsteady laminar flow about a 
rotating disc whose angular velocity w may vary with time. 
The deviation of the actual instantaneous state of the flow from 
the quasi-steady state (instantaneous steady state) is deter- 
mined. From this, a simplified criterion, ® < 0.0834 w?, is de- 
rived to define the conditions under which the flow can be con- 
sidered as quasi-steady for the purposes of shear stress and 
torque computations. Since a turbulent flow responds more 
rapidly than a laminar flow, the quasi-steady criterion found here 
should also serve for the turbulent situation. 


Introduction 


ROBLEMS involving unsteady flows comprise one 
P of the most challenging research areas in fluid 
mechanics. The mathematical difficulties inherent in 
these problems have tended to restrict analysis to 
simple flow configurations and/or to simple prescribed 
time-variations of the flow. 

Our interest here is directed to external flows. Much 
of the earlier work as reported in the texts of Goldstein, ' 
Schlichting,” and Pai* had as its goal the determination 
of the time-history of the velocity field as the body is 
accelerated from rest, either impulsively or gradually. 
In many instances, the computations were limited to 
short periods following the inception of the motion. 
Also reported are computations of the periodic response 
of the flow to imposed periodic motions of the body. 

The work of Moore‘ represented a new and in- 
genious approach to the study of the unsteady boundary 
layer. Instead of seeking the detailed time-history of 
the boundary-layer development, he focused attention 
on the state of the flow at some fixed instant of time. 
In particular, he inquired as to how much the instan- 
taneous velocity distribution in the boundary layer 
deviates from the velocity distribution which would 
exist if there were an instantaneous steady state. 
Clearly, if this deviation were small, then the steady- 
state relationships could be applied instantaneously 
and the computation of engineering results such as 
skin friction would be much simplified. (The term 
quasi-steady is used to characterize an unsteady flow 
which passes through a sequence of instantaneous 
steady states.) Moore's formulation was applied* ® 
to boundary-layer flows in which the imposed velocity 
varies in a continuously differentiable (but otherwise 
unspecified) manner with time, and in which the un- 
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steady process has been in progress for a sufficient time 
so that nearly quasi-steady conditions can exist. 

Heretofore, the study of unsteady external flows 
has been confined to laminar boundary layers over 
plane and cylindrical bodies undergoing rectilinear 
motions. In the analysis to be reported here, consider- 
ation is given to an altogether different physical situ- 
ation—namely, the flow about an unsteadily rotating 
disc located in a large body of otherwise quiescent 
fluid. Our goal is to determine the deviations of the 
actual instantaneous state of the flow from the quasi- 
steady state.t} From this, it will be possible to derive a 
criterion for distinguishing when the flow may be re- 
garded as essentially quasi-steady, thereby establishing 
when it is permissible to compute results for the un- 
steady flow by instantaneous application of steady- 
state relationships. Clearly, it would be expected that 
the extent of the deviations of the actual state of the 
flow from quasi-steady would depend both on the 
rapidity of the changes in disc speed and on the response 
characteristics of the fluid. 

A schematic representation of the system under 
study is shown in Fig. 1. The angular velocity w of 
the rotating disc may vary with time in an arbitrary, 
but continuously differentiable way. The disc radius 
is taken to be of indefinite extent. A standard set of 
cylindrical coordinates 7, ¢, z is used, and to these 
there correspond the velocity components V,, V,, V:. 
It will be supposed that the unsteady motion has been 
in progress for a sufficient time so that it is possible to 
approach to nearly quasi-steady conditions—i.e., the 
starting effects have died away. The analysis is carried 


{ During the final stages of review, there came to our attention 
a note by S. D. Nigam® which considered the fluid motion, at 
very early times, about a rotating disc which had been im- 
pulsively started from rest. This problem is not related to the 
present analysis. 


Fic. 1. Physical model and coordinates. 
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out for laminar flow of an incompressible, constant 
property fluid. 


Analysis 
Conservation Laws 


The unsteady flow field about the rotating disc is 
governed by the basic conservation principles—mass 
and momentum—and it is these which constitute the 
starting point of our study. The mathematical state- 
ment of these laws appropriate to an incompressible, 
constant property fluid may be written in cylindrical 
coordinates as follows: 


(1/r)(0/Or)(rV,) + (1/r)(OV,/d¢) + (OV./d2) = 0 
(1) 


pl(DV,/Dt) — (V,2/r)| = —(@p/dr) + 
ulv?V, — — (V,/r?)] (2) 


e[(DV,/Dt) + (V,V,/r)] = —(1/r)(0p/d¢) + 
+ (2/r°)(OV,/d¢) — (V,/r*)] (3) 


p(DV./Dt) = —(0p/dz) + (4) 
where 


D/Dt = (0/dt) + V,(0/dr) + 
(V,/r)(0/0¢) + V.(0/0z) 
V? = (0°/dr?) + (1/r)(0/dr) + 
(1/r*)(02/Og?) + (02/22?) 


Momentum conservation, Eqs. (2), (3), and (4), is 
represented by the full (unsteady) Navier-Stokes 
equations. There will be no need to invoke the simpli- 
fying assumptions of boundary-layer theory. 

These equations fully describe the physical occur- 
rences within the fluid, but to complete the statement 
of the problem, it is necessary to give the boundary 
conditions. At the disc surface, the standard viscous 
flow assumption is imposed—namely, that there is no 
fluid motion relative to the solid surface. Further, 
since the disc is impermeable to mass, the normal com- 
ponent of the velocity must vanish at the surface. 
Far from the disc, all fluid velocities must vanish aside 
from the induced axial component. A formal state- 
ment of these conditions is 


V, = 0 


0 
Ve = z=0 V (5) 
V, = 0} 


where w, the angular velocity of the disc, may vary 
with time. 

Now we turn to the task of solving this formidable 
set of partial differential equations. 


Expansion About the Quasi-Steady State 


As has already been mentioned, our goal is to in- 
vestigate the deviations of the actual instantaneous 
state of the flow from the quasi-steady state. With 
this in mind, it is natural to seek solutions for Eqs. 
(1) through (4) in the form of a series expansion about 
the quasi-steady state. 
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As a necessary prelude to this series solution, we 
briefly review some aspects of the steady-state solution 
as given by von Karman. He introduced a new 
independent variable n by the definition 


n = 2(w fy) (6a) 


where v is the kinematic viscosity. Then, for each 
velocity component and for the pressure, he assumed 
profile similarity and axial symmetry, and wrote 


F(n) = V,/rw, V,/re, 


H(n) = V./(wv)'””,  P(n) = (6b) 


With the transformation given by Eqs. (6a) and (6b), 
von Karman succeeded in reducing the partial differ- 
ential equations of the problem (Eqs. (1)—(4) with 
0/0t = 0) to ordinary differential equations. The 
solution for F, G, H, and P gave the complete steady- 
state velocity and pressure fields, as well as such engi- 
neering results as the wall shear stress and the shaft 
torque. 

Now, turning to the unsteady situation, we first ob- 
serve that the quasi-steady state is defined by Eqs. 
(6a) and (6b), where F, G, //, and P are as before, but 
w now represents the instantaneous value of the angular 
velocity. Then, to determine the deviations of the 
flow from the quasi-steady state, we write 


V,/re = f(n, Bn) = 


F(n) + BFi(n) + B2F2(n) +... (7a) 
= g(n, Bn) 

G(n) + B,G,(n) + B2G2(n) +. (7b) 
V, (wv)!/? = h(n, Bn) 

p/uw = m(n, Br) = 

P(n) + B:Pi(n) + B2P2(n) +... (7@d) 


The expansion parameters depend only on time and 
are defined by 


= &/w*,... (8) 


where = dw/dt, = d*w/dt®, etc. The first term of 
each series is associated with the quasi-steady state and 
the succeeding terms give the deviations. 

While the expansion parameters 8, can be derived in 
a mathematical manner, it is more satisfying to dis- 
play their physical meaning. Each of the 8, is pro- 
portional to the ratio of two times—namely, the time 
required for impressed changes in disc speed to diffuse 
across the velocity layer divided by a time which char- 
acterizes the rapidity with which the disc speed changes. 
Now, from one-dimensional diffusion theory, the time 
t, required for any quantity to diffuse across a layer of 
thickness 6 is proportional to 6°/y. But the thickness 
6 of the velocity layer adjacent to the rotating disc is 
proportional to (v ‘w)'/*. So, for the diffusion time, we 


find 
(9a) 


ta ~ 1/w 


On the other hand, the rapidity of the changes in disc 
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Fic. 2. The tangential velocity functions G, G, Go. 


speed is given by the following set of times: 
w/a, (w/a), (w/e), ... (9b) 


The ratio of Eq. (9a) with each of the successive mem- 
bers of Eq. (9b) yields, after some rearrangement, the 
expansion parameters B2,...of Eq. (8). If the fluid 
responds quickly to impressed changes, then the £, 
will be small and Eq. (7) properly predicts that the 
state will essentially be quasi-steady. 

Having settled on the form of the series solution, we 
now return to the governing Eqs. (1)—(4). 


Ordinary Differential Equations and Their Solution 


Turning first to mass conservation, Eq. (1), we intro- 
duce the series expansion (7a) and (7c) and note that 
0/O¢ = 0 (angular symmetry). From this, it is found 
that 


f = —(1/2)(Oh/On) (10a) 
or 
F = —(1/2)H’, F, = —(1/2)M’, 
F, = —(1/2)H2’,... (10b) 


where the primes denote differentiation with respect 
ton. Next, passing to Eqs. (2) and (3), we substitute 
the series expansions for V,, V,, and V, and group 
terms according to whether they are multiplied by 
B;, B2,.... Then, making use of Eq. (10b) to eliminate 
the f functions, we are led to the following set of ordi- 
nary differential equations: 


= HH" — (1/2)H? + a1) 
= HG'-—-H'G 
H,!” = HH," — H'H,' + + 
4GG, + + (1/2)9H" 
G,” = HG,’ = H'G, + H,G’ 
+ G + (1/2)nG’ 
= 


HH," — + H"H, + AGG: + (13) 


G2”) = HG,’ — H’'G:; + — H'G+G, 
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where now we truncate the series expansions after the 
third term. 

The boundary conditions appropriate to these equa- 
tions may be found by introducing the series (7) into 
Eqs. (5). Taking cognizance of Eq. (10b), we get 


| 


H’ = = H,! = 0277 = 0 
H,', > 0 
(14) 


Eqs. (11), (12), and (13) with the boundary condi- 
tions (14) constitute a means for determining the g 
and / functions appearing in the series (7b) and (7c). 
Further, it can be seen from Eq. (10b) that the f func- 
tions are known once the solution for the / functions 
has been obtained. 

The governing Eqs. (11) for F and G were first de- 
rived by von Karman for the steady-state problem. 
Although an approximate numerical solution was 
given by Cochran,’ it was necessary to re-solve to 
greater accuracy for the purposes of this investigation. 
Solutions have also been found for Eqs. (12) and (13). 
The computations were carried out on an IBM 653 
electronic calculator using numerical techniques out- 
lined by Gregg.’ The g, h, and f functions corre- 
sponding to these solutions are plotted in Figs. 2, 3, and 
4. With this, the deviations of the instantaneous ve- 
locity distribution from the quasi-steady state may be 
computed from Eqs. (7a), (7b), and (7c). 

From the solution, the starting values (i.e., deriv- 
atives at » = (0) essential to the forward numerical 
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FLOW ABOUT AN UNSTEADILY ROTATING DISC 


integration of Eqs. (11), (12), and (13) are found to be 
those given in Table 1. These values will be of direct 
application in the shear stress and torque computa- 
tions which are carried out later. 


TABLE 1 
Velocity Derivatives at 7 = 0 
H"™0)  He%0) Gi(0) Ge"(0) 


—1.0205 0.40967 —0.046224 —0.61592 —0.36922 0.021185 


Inspection of the differential Eqs. (11), (12), and (13) 
shows that the functions G, G;, Go, ... and H, Hy, Ho, 
. must be computed in ascending order—e.g., Ge 
cannot be computed without prior knowledge of G and 
G,. A consequence of this is that we are not easily 
able to ‘‘reach’’ terms which occur later in the series. 
Since terms beyond the third are not now available, 
we must limit our discussion to situations where the 
8, form a decreasing sequence with increasing , with 
8, less than unity.* Since the successive 8, involve 
increasing powers of w in the denominator, the require- 
ment that the sequence of 8, decrease should be satisfied 
for moderate or large angular speeds. 

The governing equations for the pressure functions 
P, P,, and P: are determined by substituting the 
series (7a)—(7d) into the final momentum conservation 
equation (4). After grouping terms, there is obtained 


P’ = H* — HH’ (15) 


P,’ = — HH,’ — 
(1/2)H — (1/2)nH’ (16) 


II 


P,' = H,"” — HH,’ — — (17) 


From these, it is seen that the pressure functions are 
computed by direct numerical integration once the 
solutions for the / functions are available. A dis- 
cussion of the boundary conditions and other aspects 
of Eqs. (15), (16), and (17) will be deferred to a later 
section devoted to pressure results. 


Torque and Shear Stress Results 


The instantaneous tangential shear stress on the disc 
surface may be calculated by application of the New- 
tonian shear formula, 


Tze = 02).=0 


After introducing the series expansion (7b) for I’, and 
taking account of the definition (6a) of n, the expression 
for r becomes 


Tee, inst. = + + 
B2G2’(0) +...) (18) 
where we have appended a second subscript to denote 


the instantaneous stress. Numerical values of G’(0), 
G,'(0) and G»’(0) are listed in Table 1. 


* This latter requirement arises because of the terms 8,2Hy, 
ete. 
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The quasi-steady shear stress 7,, ; is given by 
Tee, qs = pr(vw*)'/7G"(0) (19) 


Then, the important relationship between the instan- 
taneous and quasi-steady stresses is found by com- 
bining Eqs. (18) and (19), giving 


inst./Tzy, gs = 1 + 0.59946 (@/w*) — 
0.034396(w/w*) +... (20) 


where §, and £8, have been evaluated from Eq. (8) and 
the numerical values of Table 1 have been used. For 
positively accelerating flows (@ > (0), Eq. (20) veri- 
fies the intuitive feeling that the instantaneous stress 
should exceed the quasi-steady value. 

To overcome the shear forces exerted by the fluid, a 
torque must be supplied at the shaft of the rotating disc. 
The torque M may be computed from its definition 


M = r-(2ar dr) tz, (21) 
0 

Using the results for r,, from Eqs. (18) and (19), the 

instantaneous and quasi-steady torques are found to be 


Mins. = + 
+ B2G2'(0) + ...] (22a) 


= (22b) 
With this, the relationship between V/;,.. and .W/,, is 


Minst./Mqs = 1 + 0.59946(@/w*) — 
0.034396(@/w*) +... (23) 


which coincides with Eq. (20). 


Criterion for Quasi-Steady Conditions 


Because of the great computational simplifications 
associated with the quasi-steady assumption, it is of 
interest to derive a criterion for distinguishing those 
conditions under which the flow is effectively quasi- 
steady. To find such a criterion, we focus attention 
on the torque (or shear stress) computation and inquire 
when the quasi-steady relationship [Eq. (22b)] can 
be applied with sufficient accuracy to the unsteady 
problem. Suppose that an accuracy of 5 per cent is 


+ Shear forces on one face of the disc are considered. 
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considered sufficient. Then, using Eq. (23), it follows 
that when 


0.59946(a/w?) — 0.034396(a/w*) < 0.05 (24) 


the state can be considered quasi-steady for purposes 
of the torque computation. For moderate or large 
angular velocities, the last term of Eq. (23) will almost 
always be small relative to the preceding one. Under 
these circumstances, the quasi-steady criterion takes 
the simpler and more easily applicable form 


< 0.0834 (25) 
For example, at a rotational speed of 500 rad./sec. 


(4,775 r.p.m.), the flow will be quasi-steady for accel- 
erations up to 20,850 (rad. /sec.) /sec. 


Illustrative Example 
Suppose that the angular velocity of a rotating disc 
varies as 
w = 100(1 — e~°") 


where w is in rad./sec. and ¢ is in seconds. The torque 
(23) corresponding to this expression for w is 


Minst. _ (0.59946) 
(1 — e~*")2 


(0.034396) 
(1 3 


Att = 1/2 sec., the ratio Minst,/,, has the value 1.24. 
From either Eqs. (24) or (25), we find that the state 
is effectively quasi-steady when ¢ = 1.1 sec. 


Turbulence 


Since the response of a turbulent flow should be more 
rapid than that of a laminar flow, the criterion derived 
here for quasi-steady conditions should certainly also 
serve for the turbulent situation. 


Induced Pressure 


In considering the pressure distribution which is in- 
duced by the rotating disc, we are led to some very 
interesting and thought-provoking questions. It is 
here that a careful distinction must be made between 
physical feelings associated with flow about a finite 
disc and the somewhat harder to visualize situation of 
flow about a disc of infinite extent. 

For the finite disc, there is a finite radial discharge. 
Further, the free-stream flow which is induced to supply 
this discharge will be recruited from regions which ex- 
tend outward beyond the dise as well as from regions 
directly above the disc. So, it seems reasonable to 
expect that as we go to sufficiently great distances 
above the dise (large z), there will be felt no effect of 
the presence of the disc, regardless of whether its rota- 
tion is steady or unsteady. Hence, for such a system, 
there would be a temptation to say that the pressure 
p approaches a constant value for sufficiently large z. 
But, for the disc of infinite extent, the situation is 


JOURNAL OF THE AERO/SPACE SCIENCES—APRIL, 1960 


quite different. Here, the radial discharge is no longer 
finite, and the induced free-stream flow needed to 
supply this discharge approaches the disc with a purely 
axial velocity V.. In accordance with the principle of 
mass conservation, this axial downflow persists to in- 
definitely large heights above the disc. So, the presence 
of the infinite disc is felt at all positions in the free 
stream, no matter how large the value of z. Now, 
suppose that the rotational speed w of the disc varies. 
This will be reflected in a change—i.e., an acceleration— 
in the V,, which will be felt at all z. To satisfy the 
laws of motion, there must correspond to this acceler- 
ation a force, and this force will be a pressure gradient. 
So, it would appear that for the unsteadily rotating 
disc of infinite extent, there is a nonvanishing pressure 
gradient dp /dz for all z in the free stream. For a finite 
pressure at some finite z location, it would then appear 
that the pressure would become infinite as z > ©. 
This is theoretically reasonable since an infinite force 
is required to accelerate an infinite mass of fluid.* 

Our purpose in this discussion has been to make 
plausible the nature of the induced pressure results 
for the infinite disc and to point up the very basic man- 
ner in which they differ from those of the finite disc. 

Now, turning to the numerical results, it is first 
necessary to integrate Eqs. (15), (16), and (17). Since 
the pressure at the disc surface is unknown and the 
free-stream pressure varies with z, there appears to 
be no practfeally useful boundary condition available. 
Nevertheless, the integration can be carried out for- 
mally between the disc surface (subscript w) and some 
arbitrary point 7 (i.e., 2) in the free stream, giving 


P — P, = —0.3911 | 
P, — Py, = 0.4422 » — 0.099 >» 2 8.5 


and introducing these in the series (7d), there is ob- 
tained 


[p(n, t) — = 
—0.3911 + B,(0.4422 7 — 0.099) + 
82(0.1122 » — 0.4410), 2 8.5 (26) 


While Eq. (26) is a formally correct relationship be- 
tween the pressure at the surface and in the free stream, 
it is not a predictive expression, because neither p, nor 
p(n) are known. In the free stream, as expected on the 
basis of prior physical reasoning, there is a finite 
pressure gradient and the pressure grows unboundedly 
asn—> ©. 

In this light, the pressure results for the infinite disc 
must be regarded as primarily of theoretical interest. 
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A Three-Dimensional Linearized Analysis 
of the Forces Exerted on a Rigid Wing by a 
Shock Wave 


F. EDWARD EHLERS* ano E. M. SHOEMAKER** 
Boeing Airplane Company 


Introduction 


ly REFERENCE 1, formulas were obtained for the pres- 
sure induced on a rigid, moving plate by a plane 
acoustic shock wave striking the plate at a given angle. 
The edge of the plate was assumed to be parallel to the 
shock wave, and the analysis was therefore two-dimen- 
sional. In the present paper this restriction is re- 
moved. By suitable transformation, the case of a 
shock wave oblique to the edge of a semi-infinite plate 
can be solved in terms of the two-dimensional solution. 


The Formulation of the Boundary Value 
Problem 


As in reference 1, we shall consider that the shock 
wave is weak and that all disturbances propagate with 
a constant speed of sound c. Let us assume that a 
rigid semi-infinite plate lying in the X, Y plane of Fig. 
1 is moving in either the positive or negative X direc- 
tion subsonically with a uniform velocity U. Fig. 1 
depicts a positive velocity with the incident shock 
overtaking the trailing edge. At time ¢ = 0, the inci- 
dent shock passes through the origin 0 and the edge of 
the plate lies along the Y axis. Then at any later time 
t, the incident shock will have moved a distance ct from 
0 and the plate edge will have moved to the position 


X = Ut. The vector distance N from 0 to the shock 
front may be represented by the X, Y, Z components 
ctcosa, ctcos®, ctcosy 


where a, 8, and y are the angles the normal to the shock 
front makes with the X, Y, and Z axes, respectively. 
If R = (X, Y, Z) is a vector representing an arbitrary 
point on the shock plane, then it must satisfy 


(R-—N)-N=0 
In terms of the coordinates this becomes 
X cosa + Ycos8 + Zcos y = ct (1) 


To find the intersection of the incident shock with 
the edge of the moving plate, point 0’ of Fig. 1, we set 
X = Utand Z = Oin Eq. (1). The resulting solution, 
denoted by Y, becomes 
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Fie. 1. 


Y,/ct = [1 — Moy cos a)/cos (2) 


where My) = U/c is the Mach number of the half plane. 

The locus of these intersection points in the XY, Y 
plane describes a straight line which makes an angle 
6, with the X axis shown in Fig. 1 and given by 


0, = tan" Y./Ut = 
tan~! [((1 — My cos a)/My cos B] (3) 


This line is also the axis of the cone formed by the 
envelope of the disturbances produced by the shock 
front intersecting the edge of the plate. Since all 
disturbances propagate with the velocity c, the cone 
is tangent to the shock front along the line passing 
through the point of intersection (Ut, VY.) and the 


intersection of N with the shock plane. The cone 
generating angle is 
6. = sin! [ct, V (Ut)? + (Y,)?] | 
(4 
= sin- (cos 


The regions outside of the cone are regions of con- 
stant state with known pressures. The surface of the 
cone is divided into sections of constant pressure by 
rays through the vertex; a typical cross section in a 
plane normal to the cone axis is illustrated in Fig. 2. 
On the section of the cone surface represented by the 
larger arc FH the pressure has the same value as the 
pressure behind the shock before it strikes the plate. 


The wave pattern for a plane acoustic shock wave strik- 
ing a half plane traveling subsonically. 
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REG (0) 
7 REG (1) 
\ 
H 
Fic. 2. A typical cross section in a plane normal to the cone 
axis 00’ of Fig. 1 for a shock front overtaking a plate traveling 


subsonically. 


When the pressure is normalized to the magnitude of 
the jump in pressure across the incident shock, the 
value of the pressure P on FH is unity. Since the 
disturbance FB+ is propagating into the region un- 
disturbed by the plane shock, the value of the normal- 
ized pressure on FB+ is zero. Similarly, the disturb- 
ance on are /7B— is propagating into the region behind 
the reflected shock, and the value of the pressure on the 
part of the cone represented by //B— is equal to 2. 
There remains the determination of the pressure field 
interior to the cone, in particular, on the plate repre- 
senting the wing. 

The pressure field inside the cone is constant on rays 
through the vertex of the cone. The pressure is then a 
function of two variables represented by the conical 
coordinates with the vertex of the cone as origin. It is 
convenient to take advantage of the conical nature of 
the field whereby the three-dimensional field equations 
can be reduced to the two-dimensional case of refer- 
ence 1. 


Derivation of the Differential Equation in 
Conical Coordinates 


If we move with the reference frame &, 7, ¢ having 
its origin at the cone vertex then the geometry and the 
pressure field are stationary in time and depend only 
upon (é, , ¢). We shall find first the coordinate trans- 
formation which transforms the X, Y, Z axes to &, 7, 
¢ axes with origin on the cone vertex and the 7 axis 
oriented along the axis of the cone. The Z axis remains 
unchanged and the &, 7 coordinates are defined by 


& = (X — Ut) sin®@, — (Y — Y,) cos Ba | 
n = (X — Ut) cos, + (Y — Y,) sin 6, } (5) 
With Y, = Ut tan @,, the equations for the coordinates, 
£, n, ¢, simplify to 
= Xsin@, — Y cosh | 
n = X cos0@, + Y sin @, — Ut sec 0, (6) 
| 


The pressure field satisfies the wave equation 


Pxx + Pyy + Pez — P.u/c? = 0 (7) 
Transforming to the variables é, n, ¢, Eq. (7) becomes 
Pe + Pry — (Mo? sec? 0, — 1)P,, = 0 (8) 
From Eq. (4) we find that 
M,? sec? 6, — 1 = cot? 6, 
so that Eq. (8) reduces to 
Pe + — P,,/tan? 0, = 0 (9) 


This is the wave equation in two dimensions with 
n tan 6, replacing ct. 
Introducing the conical transformation 


x = —&/ntan6,, y = —f£/n tan 6, (10) 


and the Busemann transformation 


mn =r/(1+ (1 wherer = 
Vx? + y?, (11) 
= (y/x) 


reduces Eq. (9) to Laplace’s equation. Under these 
transformations the cone is transformed to the unit 
circle in both the x, y plane and the (r, @) plane. We 
have therefore reduced the problem to the two-dimen- 
sional one discussed in reference 1. 


The Parameters Associated With the 
Two-Dimensional Solution 


There remains the task of finding the values of 
and @, associated with the two-dimensional problem of 
reference 1 in terms of the variables Mo, the actual 
plate Mach number, and the direction cosines of the 
normal to the plane shock front. The edge of the plate 
is given by ¥ = Ut. In the (é, n) plane this is 


—(Y — Y,) cos @ 
= (Y — Y,) sin 


ll 


In the x, y plane this becomes 
tan 0, = cot 6, cot | 
f (12) 
Since in the x, y plane the cone is a circular cylinder of 
unit radius, the plate Mach number is J = z. There- 
fore, substituting for 6, from Eq. (4) in Eq. (12) yields 


M = V (Mo? — cos? @,)/(1 — cos? @) (13) 


The angle 6; is determined by the intersection of a 
plane normal to the cone axis with the &, » plane and 
with the plane determined by the cone axis and the 
normal vector N. The vector normal to the plane 
defined by the cone axis and the normal vector N is 
C X N where C is a unit vector along the n axis, The 
components of the vector N in the &, 7, ¢ system are 
found by putting 


X =ctcosa, Y=ctcos®, Z=ctcosy 


in Eq. (6). The & », ¢ components of N then are pro- 
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portional to 
cos a sin 6, — cos 8 cos@ = & 


cos a cos 8, + cos B sin 0g — Mo sec @. = m 


cos y = %1 
and C is given by 
(0, 1, 0) 
Hence C X N = (i, 0, —£&) and the angle 6, is given by 
6, = — tan! (&/f1) (14) 


or 


6, = 2/2 — [(cos sin — cos B cos @,)/cos 
(15) 


When @, is eliminated by Eq. (3), we get 
6, = 2/2 — tan! [(cos a — Mo sin® y) + 


V1 + — 2M, cos a — “Me? cos? (cos y)] (16) 


General Properties of the Pressure Distribution 


From the preceding analysis and the properties of 
the two-dimensional solution, we can infer some of the 
properties of the pressure distribution on a wing. 
The shock moves along the trailing edge of the wing 
with a velocity derived from Eq. (2) of 


(1 — My cos a)c/cos B (17) 


From this we see that the shock wave does not move 
over the wing if 


Moy = 1/cos a 


Note that this is a supersonic Mach number. For 
all Mach numbers less than 1/cos a, the shock wave 
moves along the edge of the wing and eventually 
envelops the entire wing. The resulting pressure dis- 
tribution introduces a twisting moment on the wing 
when 8 # 7/2. 


PLATE MOTION 


u— 


ox 


Fic. 3. The net pressure distribution on a semi-infinite plate 
traveling subsonically at a Mach number of 0.4 and overtaken 
by an acoustic shock front having direction cosines a = 45°, 
8 = 60°, and y = 60° (see Fig. 1). Isometric drawing. 
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Fic. 4. A typical cross section in a plate normal to the cone 
axis for a plate traveling supersonically toward the shock front. 


From Eq. (13), the cone envelope of disturbances is 
tangent to the edge of the plane when ./ is equal to 
unity in Eq. (13). The actual Mach number My of 
the plate then also is unity. The pressure field then 
is quite simple. On top of the plate the pressure is 
identically zero. On the lower side the pressure is 2 
behind the shock wave and zero elsewhere. 

As a numerical example, we consider a flat plate 
traveling in the positive X direction at a Mach number 
of 0.4. The normal to the shock front is assumed to 
make the angles 


a= 45°, B= 60°, = 60° 


with the X, Y, Z axes, respectively. The angle 6, be- 
tween the X axis and the axis of the cone is calculated 
from Eq. (3) and found to be 74°25’. The cone gener- 
ating angle obtained from Eq. (4) is 42°11’. To cal- 
culate the actual pressures from Eqs. (3.12) to (3.17) of 
reference 1, we must evaluate @, from Eq. (15) and 7 
from Eq. (13). These values are 


6, = 42°27’ and M = 0.306 


The resulting distribution of the difference in pressure 
between the upper and lower surfaces of the plate is 
shown plotted along lines perpendicular to the cone 
axis in Fig. 3. The twisting moment at each cross sec- 
tion produced by the pressure distribution varying 
along the plate span is well illustrated by the drawing. 


The Shock Wave Striking the Leading Edge of a 
Supersonically Moving Plate 


The three-dimensional problem of the shock wave 
striking the leading edge of a supersonically traveling 
plate also can be reduced to the corresponding two- 
dimensional problem discussed in Section (3) of refer- 
ence 1. The cone of the envelope of disturbances from 
the intersection of the edge of the plate and the shock 
wave is cut completely by the plate. This is illustrated 
by the cross section in a plane normal to the cone axis 
in Fig. 4. 

To find the angle uw needed in the two-dimensional 
solutions, we require the coordinate ¢ of the wing lead- 
ing edge, which is given in Eq. (12). The axis angle 6, 
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is given by Eq. (3) with —M replacing Mo. Thus, 
6, = tan~! [(1 + My cos a)/(— Mo cos B)| (18) 


For a and 8 < 7/2, 6, is an angle greater than 7/2 and 
the leading-edge coordinate # is therefore negative. 
This coordinate may be expressed in terms of 6, and 


My by using the relation following Eq. (8). Thus, 
= = —cos? 0,/sin 0,.V My? — cos 9% (19) 


Since the cone becomes a unit circle by the conical 
transformation, the Mach number J/ to be used in 
defining the angle ~« = cos~! (1/1) is given by the 
magnitude of x coordinate of the leading edge or \z | of 
Eq. (19). 

We now must find the magnitude of the pressure on 
axis B+’F’ and B-'H’ in Fig. 4. This is a region of 
steady-state uniform supersonic flow, and we can make 
use of linearized theory. Let u be the magnitude of 


0) = + 0, + (1/M) — 2 


and 


the perturbation velocity behind the undisturbed plane 
shock, then the angle of attack ap of the flow in the 
leading-edge region | B of the plate is 


ay = tan-! [(u cos y)/(U + u cos a)] | 


(20) 


Following the same analysis as in Section (3) of refer- 
ence 1, we obtain the following values of the normalized 
pressures in regions (3) and (4), respectively, of Fig. 4: 


ay  (u cos y)/U 


P3 1 — Mo cos Me? — 1 (21) 


Ps = 1+ Mocos Me? — 1 (22) 


Comparing Eqs. (21) and (22) with Eqs. (3.6) and (3.7) 
of reference 1, we see from Eqs. (3.12) and (3.14) of 
reference 1 that the pressure inside the unit circle on 
the top of the plate is given by 


P=P,+ (23) 


where 


(1 + Mx,) sin + (cos — x1) Vie 


(24) 
(1 + Mx,)(cos — x) — VM? — 1 sin A; 


P(x, 0) = (1 — My cos y/V Mo? — 1){2 — 1/(1 + Mx,)] — cos-! (25) 


The Mach number J/ to be used in Eqs. (24) and (25) 
is given by |z| of Eq. (19) and the angle @, by Eq. (16) 
with @, defined by Eq. (18). The coordinate x; is re- 
lated to the conical coordinate x by Busemann’s trans- 
formation 


= x/1t+ x?) 


Conclusion 


The pressure distribution on a moving flat plate in- 
duced by an acoustic shock front striking the edge of 
the plate obliquely has been found in terms of the two- 
dimensional solution of reference 1. The parameters 
required by the two-dimensional solutions are given 
as functions of the Mach number of the plate and of 
the direction cosines of the normal to the shock front. 
These direction cosines are referred to an orthogonal 
coordinate system oriented with one edge of the plate, 


with the direction of motion of the plate, and with the 
normal to the plate surface. The pressure at any in- 
stant is constant along rays through the point of inter- 
section of the shock front and the plate edge. How- 
ever, the pressure has not been found in the region of 
the corner of the plane contained in the sphere of the 
disturbances from the shock wave just touching the 
corner of the plate. The pressure here would be some- 
what less than that predicted by the conical solution 
because of the pressure relief provided by the corner. 
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The Effect of Controlled Three-Dimensional 
Roughness on Boundary-Layer Transition 
at Supersonic Speeds' 


E. R. van DRIEST* ano W. D. McCAULEY** 


North American Aviation, Inc. 


Summary 


Experiments were performed in the 12-in. supersonic wind 
tunnel of the Jet Propulsion Laboratory of the California Insti- 
tute of Technology to investigate the effect of three-dimensional 
roughness elements (spheres) on boundary-layer transition on a 
10° (apex angle) cone with zero heat transfer. Data were ob- 
tained at local Mach numbers of 1.90, 2.71, and 3.67 by varying 
trip size, position, spacing, and Reynolds number per inch. 
The results indicate that (1) transition from laminar to turbulent 
flow induced by three-dimensional roughness elements begins 
when the double row of spiral vortices trailing each element con- 
taminates and breaks down the surrounding field of vorticity, (2) 
transition appears rather suddenly, becoming more violent with 
increasing roughness height relative to the boundary-layer thick- 
ness, (3) after the breakdown of the vorticity field, the strength 
of the spiral vortices may still persist in the sublayer of the en- 
suing turbulent flow, (4) lateral spacing of roughness elements has 
little effect upon the initial breakdown (contamination) of the 
laminar flow, and (5) the trip Reynolds number (uk/v)s5, where u 
and v are the velocity and kinematic viscosity at the outer edge of 
the boundary layer and k is roughness height, such that transition 
occurs at the roughness position, varies as the position Reynolds 
number to the one-fourth power, viz. (ux,/v)g!/4, where x, is 
trip position. 


Symbols 
x = distance along ray from point of cone 
k = roughness height 
Xo = distance to transition for a smooth cone 
xr = distance to transition for a tripped boundary layer 
Xk = distance to trip 
u = velocity in x direction 
p = fluid density 
m = coefficient of viscosity 
v = kinematic viscosity 
T = shear stress at wall 
Reo = Reynolds number of transition for smooth cone 
Rer = Reynolds number of transition for tripped boundary 
layers 
Re; = Reynolds number of trip 
Re,, = Reynolds number of trip position 
Ms = Mach number at outer edge of boundary layer 
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T,../Ts = wall to free-stream temperature ratio of boundary 
layer 

r = recovery factor of boundary layer 

6 = subscript denoting conditions at outer edge of 


boundary layer 


Introduction 


™ PHENOMENON Of transition from a laminar to a 
turbulent boundary layer has been the subject of 
intense research in recent years to learn the basic mecha- 
nism of the transition process. A knowledge of factors 
controlling transition is, of course, directly applicable 
to the design of aircraft such as high-altitude glide 
vehicles and hypersonic re-entry bodies, particularly be- 
cause of the role transition plays in aerodynamic heat- 
ing. 

The present paper describes experimental studies of 
transition at supersonic speeds without heat transfer, 
using different-sized three-dimensional surface rough- 
nesses (spheres) to disturb the flow. The experiments 
were conducted on a 10° (apex angle) cone in the 12-in. 
supersonic wind tunnel of the Jet Propulsion Labora- 
tory of the California Institute of Technology. The 
roughnesses were located at the 5-, 7-, and 9-in. sta- 
tions, and the Reynolds number per unit length of the 
tunnel was varied to extend the range of effectiveness 
of the roughnesses. Transition was observed by three 
complementary methods—namely, (1) anamorphic 
schlieren, which magnified the image of the boundary 
layer twenty times normal to the flow, (2) surface tem- 
perature recovery measurements, and (3) surface flow 
pattern indicators. The local Mach numbers at the 
outer edge of the boundary layer were 1.90, 2.71, and 
3.67. 

The present work is an extension of previous studies 
utilizing single two-dimensional roughness elements 
(wires) to promote transition.! 


Apparatus 


Wind Tunnel 


The experiments were conducted in the 12-in. super- 
sonic wind tunnel of the Jet Propulsion Laboratory 
(CIT). The tunnel is of the closed circuit, variable 
pressure, continuous operation type. Since the Mach 
number could be varied from approximately 1 to 4 
through the use of a flexible-wall nozzle, the present 
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tests were run at free-stream Mach numbers of 1.97, 
2.81, and 3.84. These Mach numbers correspond to 
local Mach numbers of 1.90, 2.71, and 3.67, respec- 
tively. During the tests the supply temperature varied 
from approximately 90°F. to 120°F. with a dew point 
normally below —10°F. The Reynolds number per 
inch at local Mach number 1.90 could be varied from 
0.3 X 108 to 0.5 & 108; at Mach number 2.71 from 0.2 
xX 10° to 0.7 & 10°; and at Mach number 3.84 from 
0.2 X 10° to 0.56 X 108. The normal operating turbu- 
lence level of the tunnel settling chamber was about 
0.3 per cent. 


Test Body 


The test body used in this investigation was the 
same as that of reference 1, being a 20-in. 10° (apex 
angle) cone with a 321 stainless steel shell 1/8-in. thick 
(Fig. 1). The shell was instrumented with 35 stainless 
steel-constantan thermocouples located along a ray on 
the outer surface of the cone and electroplated with a 
0.003-in. coat of electrolus nickel. For these tests the 
tip radius was less than 0.0005 in., and the root mean 
square roughness was estimated to be less than 10 
microinches. 


Temperature Instrumentation 


The thermocouples were formed using the cone as a 
common stainless steel leg and constantan wires silver 
soldered to the cone as the other legs. The constantan 
leads and a single stainless steel wire from the cone were 
passed through a multichannel Ledex automatic selec- 
tor switch to reduce the number of leads to a single pair 
of stainless steel and constantan leads. The single con- 
stantan wire then went to an ice-bath junction from 
which it emerged as a stainless steel wire. This latter 
wire and the common stainless steel lead from the cone 
were then led into a Brown Potentiometer recorder 
which was equipped with a Gianini digitizer. The digi- 
tizer data were fed into a punch-tape machine, thence 
to a computer. The reduced data came from the com- 
puter on a tape from which the data were automatically 
typed and plotted. 


Magnified Schlieren System 


A special optical system,! which magnified the 
schlieren image of the boundary layer by a factor of 20 
normal to the flow, was used to observe the boundary 


layer. 


Discussion of the Experiments 
(a) Transition on a Smooth Cone With Zero Heat Transfer 
At a given Mach number, transition occurs on a 
smooth cone as a result of the dynamics of the flow and 
external (free-stream) disturbances. Thus, the Rey- 
nolds number of smooth-wall transition is described 
functionally by 


= f(M;, dimensionless disturbance numbers) 


(1) 
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in which p;, “5, us, and M; are the density, velocity, 
viscosity, and Mach number, respectively, of the local 
flow at the outer edge of the boundary layer, and x is 
the distance from the point of the cone to the end of 
transition. For constant Mach number and similar 
disturbances, Eq. (1) then becomes 

PssXo/us = Constant (2) 
Hence, if dynamic similarity prevails, the distance to 
transition x) should vary inversely with p;u5;/u;. The 
factor psu5/us is called the ““Reynolds number per unit 
length” or simply “‘unit Reynolds number.” 

In the present experiments, the unit Reynolds num- 
ber is varied arbitrarily (by changing the density p; 
through change in pressure level) to give a large range 
of roughness Reynolds number over which the effect of 
roughness on transition can be studied. 

Fig. 2 shows the distribution of the wall-to-local-free- 
stream temperature ratio 7,,/7; along the cone for dif- 
ferent Reynolds numbers per unit length and local 
Mach number 2.71. The corresponding magnified- 
schlieren photograph is attached above each set of 
data. It is immediately noted that (1) the schlieren 
shows scuffing up of the laminar boundary layer just 
before the peak of the temperature distribution (the 
initial rapid thickening of the boundary layer is defined 
in this paper as ‘“‘transition’’), (2) the distance to 
transition x9 increases as the unit Reynolds number de- 
creases in accordance with Eq. (2), (3) the slope of 
the temperature rise decreases with increase in transi- 
tion distance, and (4) the temperature overshoots at 
transition. 

Since the recovery factor r is a better parameter than 
T../T; for studying the behavior of transition with 
zero heat transfer, because the recovery factor is prac- 
tically independent of Mach number whereas 7,,/7; 
is not, the data of Fig. 2 are replotted in Fig. 3 in terms 
of recovery factor. Of significant importance is the 
decrease in slope of the recovery factor (temperature) 
rise with increasing distance to transition, which indi- 
cates an increasingly longer local breakdown of the 
laminar flow. 

It is now instructive to change the scale of the ab- 
scissa of Fig. 3 from length x to local Reynolds number 
pssX/us. This is done in Fig. 4. That the transition 
regions now fall close together indicates that the flow 
phenomenon remains in effect similar when p,u5/us is 
changed. This demonstrates that on a Reynolds num- 
ber scale the breakdown of the laminar flow is ac- 
complished in the same relative manner. 

The overshoot in recovery temperature at transition 
is undoubtedly due to the sudden collapse of the vorti- 
city field before a statistical turbulent regime predomi- 
nates. It indicates that transition is a more violent 
phenomenon than previously expected, particularly in a 
low turbulence wind tunnel. 


(b) Transition With a Flat Band Mounted on 
the Smooth Cone 


In order that different three-dimensional roughnesses 
might be readily placed on the model during testing, 
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Fic. 2. Magnified schlieren pictures and temperature distributions for a smooth adiabatic wall. M/, = 2.71. 
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two-dimensional bands approximately 0.001 in. thick 
and 0.25 in. wide were constructed to mount the various- 
sized spheres on. This method of mounting rough- 
nesses was chosen, because previous tests! had shown 
that two-dimensional wires of such height would not 
affect smooth-wall transition for the Mach numbers of 
the present tests. Comparator pictures of typical ball- 
band combinations mounted on the model are pre- 
sented in Fig. 5, where it is seen that the fillets between 
the spheres and bands were small and the bands fitted 
the model snugly. 

A series of tests was first carried out in order to ensure 
that the bands alone would not alter the smooth-body 
transition position. It was found from schlieren photos, 
surface temperature measurements, and azo-benzene 
patterns that the band alone did not appreciably affect 
the smooth-wall data for MM; = 2.71 and 3.67. It was 
therefore concluded that the bands were thin enough to 
mount sufficiently large three-dimensional roughnesses 
on them for transition studies at Mach numbers 2.71 
and 3.67. For precision, direct mounting of spheres on 
the cone were necessary at 7; = 1.90. 


(c) Transition From Three-Dimensional 
(Spherical) Roughnesses 

Fig. 6 shows typical transition data as indicated by 
magnified schlieren, azo-benzene, and temperature dis- 
tribution for 0.010-in. spheres mounted at 1/16-in. 
spacing around a band at the 5-in. station for flow at 
Mach number 2.71 with zero heat transfer. The 
smooth-cone temperature results of Fig. 2 are also 
plotted in the figure. 

As transition moves back on the model, each method 
of transition observation has advantages in a different 
region. Although all of the methods agree reasonably 
well, the regions of best use are: (1) when transition is 
close to the trip, the azo-benzene pattern gives a sharp 
determination of the transition; (2) as transition 
moves further back on the model, the azo-benzene pat- 
terns appear to wash out due to nonuniform spraying 
on top and bottom, and the magnified schlieren picture 
yields a more definite delineation; and (3) when 
transition is far back on the model, the magnified 
schlieren also becomes blurred, and the temperature 
peak becomes more reliable for determination of the 
transition point. Finally, for quantitative purposes, 
transition data were obtained after checking each 
method against the others. 

The temperature distributions of Fig. 6 warrant some 
discussion. At the Reynolds numbers per inch such 
that the roughness pulls transition forward almost to the 
trip, it can be noticed that the temperature rises even 
ahead of the roughness elements. Since the laminar 
temperature distribution for the smooth-wall transition 
did not begin to rise in the region ahead of the roughness 
position, it appears proper to conclude that the tem- 
perature rise ahead of the roughness may be attributed 
to longitudinal heat conduction in the model shell. 
As the Reynolds number per inch is decreased, the slope 
of the temperature distribution between the trip and 
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Fic. 7. Comparison of transition due to three different roughness 
diameters. M,; = 1.90. 


transition becomes more and more shallow. Part of 
the decrease in slope may be attributed to the increase 
in distance to transition. However, the remainder of 
the temperature flattening is undoubtedly due to the 
vortex paths induced by the roughness spheres. 

In order to obtain the effect of Mach number, further 
data were taken at local Mach numbers 1.90 and 3.67. 

Fig. 7 for 17, = 1.90 is very instructive, because it 
illustrates in a comparative way the control of transi- 
tion by different sized trips as the thickness of the 
boundary layer—i.e., unit Reynolds number—is 
changed. 

Fig. 8 is typical of 17, = 3.67. It is interesting for 
two reasons: it shows (1) that the eroded path behind 
each roughness element is really a double row of longi- 
tudinal vortices with a stagnation line between them, 
and (2) that when initial contamination sets in—i.e., the 
laminar flow breaks down—the process is immediate 
and definite, forming fan-shaped patterns in the azo- 
benzene.t 


t+ It has been brought to the authors’ attention that Gregory 
and Walker of the NPL (R&M 2779) made similar observations 
at subsonic speeds of the roughness induced surface pattern be- 
hind cylindrical elements and drew the same conclusion concern- 
ing the spiral vortices trailing the roughness elements. 
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The direction of the vortices trailing each roughness 
element is shown in Fig. 9. This follows from the con- 
cept that the vorticity of the oncoming boundary layer 
“‘wraps itself around” the roughness element. Note that 
scouring would also occur in front of as well as around 
the sides of the spherical element. 


(d) Effect of Lateral Roughness Spacing on Transition 


In many of the tests, the lateral spacing of the rough- 
ness elements was varied from 1/16 to 1/8 to 3/16 in. 
in order to determine whether the spacing had any in- 
fluence on transition. This also allowed closer scrutiny 
of the azo patterns. 

From the tests with the wider spacing, it was found, 
as seen in Fig. 8, that breakdown of the flow is asso- 
ciated with a fan-shaped zone of local turbulence, whose 
angle of lateral spread is about 5° at W/; = 1.90, 6° at 
M, = 2.71, and 7° at WM; = 3.67. Such spreading is 
called transverse contamination after Charters.” 

As far as spacing is concerned, the tests showed that 
the position of initial contamination is independent of 
the spacing of the roughness elements, provided the ele- 
ments are not too close together—i.e., provided the 
trailing rows of double spiral vortices do not interfere 
with each other. In these experiments, the 1 16-in. 
spacing was still sufficiently large to cause no inter- 
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Double row of spiral vortices behind each sphere in 
boundary-layer flow. 
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ference with the neighboring vortex patterns uutil the 
fan-shaped contamination areas met. 


(e) Effect of Longitudinal Position on Transition 


In order to study the effect of longitudinal position of 
the roughness elements on transition, tests were also 
conducted with elements at the 7- and 9-in. station for 
M; 


(f) Miscellaneous Information From the Observations 
(1) Explosive Nature of Transition—The distinct 
thickening of the boundary-layer schlieren photos and 
the overshoot of the temperature profiles of Fig. 2 indi- 
cate that the process of transition, certainly in the later 
stage, is more violent (in fact almost explosive) than 
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Fic. 8. Comparison of transition due to controlled roughness by magnified schlieren, azo-benzene patterns, and temperature distributions 
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TRIP CONT. AMINATION 


(a) LOCAL SCOURING AT TRANSITION (AZO-BENZENE ) 
Mg=567, Xx =5", k=O0I5", BALL SPACING = 1/32" 


TRI ad CONTAMINATION 


(b) PERSISTENCE OF VORTICES IN TURBULENT SUBLAYER (OIL FLOW; 
AZO-BENZENE INSET) 


Mg= 271, x, =5",k=0.015", BALL SPACING = 1/16" 


TRIP CONTAMINATION 


(c) INDUCTION OF LATERAL VORTICES IN CONTAMINATION REGION 
(OIL FLOW) Mg=2.71, x, #7",k* 0010", BALL SPACING = 3/16" 


Fic. 10. Surface flow visualization of the roughness induced transition process. 
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Fic. 11. Magnified schlieren pictures and ——— distribu- 
tions showing laminar sublayer. M,; = 
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Fic. 13. Effect of Reynolds number per inch on transition 
induced by three-dimensional roughness at 5-in. station. M; = 
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Fic. 14. Effect of Reynolds number per inch on transition 
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had been anticipated. This phenomenon is also shown 
very clearly in the scoured-out portion of the azo-ben- 
zene pattern of Fig. 10(a) obtained at 1/, = 3.67 with 
spherical trips of 0.015-in. diameter spaced at 1/32 in. 


(2) Persistence of Longitudinal Vortices in Sublayer of 


Turbulent Flow—Close observation of Fig. 10(a) shows 
another interesting phenomenon—viz., that the spiral 
vortices trailing the roughnesses still persist in the sub- 
layer of the ensuing turbulent flow even after they have 
completely entangled the vorticity in the laminar layer 
and caused transition. In order to verify this sublayer 
pattern, separate experiments were performed using an 
oil-flow visualization technique in which the model is 
covered with a thin coat of SAE-30 mineral oil and ob- 
served under ultraviolet light during the test. Fig. 
10(b) is a composite of azo-benzene and oil patterns 
showing the surface patterns through transition region. 

(3) Development of Contamination Region—Fig. 10(c) 
illustrates an oil-flow pattern for the entire region behind 
a ring of 0.010-in. spheres spaced 3/16 in. apart with 
M; = 2.71. It is interesting to observe that in the 
region of contamination the double row of vortices give 
rise to further lateral vortices which can be seen to dis- 
sipate downstream. Note that the sublayer pattern in 
Fig. 10(c) dissipates quicker than in Fig. 10(b), be- 
cause the trips in Fig. 10(c) are smaller and are stationed 
further downstream. 


(4) Observation of Viscous Sublayer of Turbulent 
Flow—In the course of enlargement of some schlieren 
photos for flow at low unit Reynolds number, it was 
noticed that the outline of the boundary layer began 
to thin out rather than thicken when transition was 
supposed to set in as indicated by the temperature 
distribution. On second thought, it was realized that 
the thinned-out portion of the schlieren photo was un- 
doubtedly the viscous sublayer of the turbulent 
boundary layer, and that what had happened was that, 
in decreasing the unit Reynolds number, the density of 
the flow was becoming so low that the schlieren could 
no longer discern changes in the turbulent part of the 
layer, but only in the sublayer, where in fact the great- 
est density gradient occurs. Fig. 11 shows the sublayer 
shifted forward as the temperature peak (transition) 
is made to move forward by roughness elements at the 
5-in. station. 


(5) Microsecond Shadowgraphs—Microsecond shad- 
owgraphs of the boundary layer with transition induced 
by roughness were taken and compared with the steady- 
state transition position, as determined from the long- 
period magnified schlieren pictures. It was observed 
that the steady-state location of transition agreed well 
with the instantaneous location. It was next noted that 
transition occurred symmetrically around the cone. A 
further observation was that the flow was essentially 
smooth before transition—i.e., large local turbulent 
bursts of the type observed by Emmons® were not de- 
veloped until possibly late in the transition process. 
These observations all attest to the smoothness of the 
free-stream flow in the supersonic wind tunnel, which 


smoothness of flow is absolutely essential to the success- 
ful conduct of transition studies. 


Correlation of the Data 


Measurements of transition position x; obtained in 
the manner described above are first plotted versus unit 
Reynolds number in Figs. 12, 13, and 14 for the three 
local Mach numbers 1.90, 2.71, and 3.67, respectively. 

These plots show immediately how well Eq. (2) 
(which represents a hyperbola) is realized for the 
smooth-wall case. In fact, at JJ; = 1.90, the smooth- 
wall data agree exactly with Eq. (2), indicating effec- 
tive dynamic similarity as far as free-stream dis- 
turbances are concerned; however, the data deviate 
increasingly with increased Mach number. 

Next of interest is the rapid movement of transi- 
tion as the “effective roughness size’’ is approached. 
“Effective roughness size’’ is defined in this report as 
that size required just to bring transition forward to the 
trip. It is further noted that the distance between the 
trip and the transition position when the trip is effective 
increases with Mach number. 

The data of Figs. 12-14 are replotted in Figs. 15-17. 
Extending Eq. (1) to include a row of single-element 
roughnesses of height k located at position x, on the 
cone, one obtains by dimensional analysis for induced 
transition at position x; 


= us, Ms, 
dimensionless free-stream disturbance numbers) (3) 


Since the intensity of disturbances in the free stream 
establishes the smooth-wall transition Reynolds num- 
ber corresponding to position x» by Eq. (1), Eq. (3) may 
be written as 


bs = 
Ms, ur, Ms) (A) 


Furthermore, if remains appreciably con- 
stant during experiments with roughness, and this is 
seen above to be approximately the case, at least at the 
lower Mach numbers; then, for each Mach number, Eq. 
(4) would become 


P us = Ms, us) (5) 


A small effect of psu3xo/us; may however be taken into 
account by a linear proportion, and therefore Eq. (4) 
would be 


Rer/ Rey = f(Rez,, Rex) (6) 
for some average value of Reo, where 
Rer = psusXr/us, Reo = psttsXo/ us, 
Rez, = psltsXi/us, Rex = pattshk/us 


For this reason, the ordinates of Figs. 15-17 are normal- 
ized by dividing by Reo, and in each figure an average 
value of Reo is given. 

The general effect of trip position Reynolds number 
Re,, is not explicitly indicated in Figs. 15-17. However, 
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it can be accounted for in the special case when the 
roughness Reynolds number Ke, has such a value that 
transition is brought forward just to the roughness ele- 
ment—i.e., Rep = Re,,. For such a condition, Eq. (6) 
reduces to 


Re,,/Rey = f( Rex) (7) 


which then gives the minimum height of spherical 
roughness element required to induce transition at the 
trip position. The transition Reynoods number is also 
then a minimum. Such a function which appears to fit 
the minimum values of Rer/ Rey is plotted as the dashed 
curve in each of Figs. 15-17. 

The form of Eq. (7) may be deduced as follows: It 
was pointed out in reference 4 that the size k of three- 
dimensional roughness element required to break up the 
viscous sublayer of a turbulent flow along a smooth wall 
and produce “rough’’ conditions was dictated by the 
expression evs ‘p k/w = 60, where 7 is the local shear 
stress. Thus, k varies inversely as the square root of 
the shear stress—i.e., 


k~1/Vr (8) 


or, assuming a 1/7 power velocity profile, which implies 
a 1/5-power skin-friction variation 


k~ (9) 


where z is the distance from transition to the point in 
question in the turbulent boundary layer. 

If it is assumed that Eq. (8) can be applied to the 
laminar flow region preceding transition, then the 
three-dimensional roughness height k required to cause 
immediate transition should vary as 


Re, ~ (10) 


which is the form of the dashed curves in Fig. 15-17. 
This thickness k might be called the ‘‘stability sub- 
layer” of the laminar flow. 

If, furthermore, Eq. (7) is general, then it should 
hold when the roughnesses are placed at different sta- 
tions. To check this effect, roughnesses were also 
placed at the 7- and 9-in. stations for M; = 2.71. The 
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Fic. 22. Minimum trip Reynolds number at which transition 
Reynolds number is 1,000,000 on a 10° cone as a function of Mach 
number. 
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data are shown in Figs. 18 and 19. Apparently Eq. 
(10) still holds, as evidenced by the dashed curve in 
Fig. 19. 

The inner stability regions of the laminar and turbu- 
lent boundary layers are illustrated schematically in 
Fig. 20. 

When the dashed curves of Figs. 15-17 are trans- 
posed to a single sheet, Fig. 21 results. The incom- 
pressible-flow curve of Fig. 21 was estimated by as- 
suming pu;6/u = 540 as the beginning of the stability 
sublayer of the laminar boundary layer. 

A cross section through Fig. 21 at Rey = 10° yields 
Fig. 22, which tells what size spherical roughness is re- 
quired just to hold transition at the roughness element 
with a value of Rey = Re, = 10°. Curves for other 
Reynolds numbers can be similarly plotted. The prac- 
tical use of this information for inducing transition 
in models testing is immediately apparent. 

Eq. (10) may be written in terms of k/6* where 6* is 
the boundary-layer displacement thickness at the trip. 
Thus, 


Ren ~ (k/8*)~4 (11) 


This equation is of interest because of.the method used 
by Dryden? in correlating the transition data of Tani, 
Hama, and Mituisi® for single wire trips on a flat plate 
at low speeds. While Dryden’s analysis indicated that 
k /6* was sufficient to correlate all transition data from 
smooth wall to tripped transition, Eq. (11) above is 
held valid at least for the case where transition just 
meets the trip. 

Fig. 23 is finally inserted to show the relative effects 
of two- and three-dimensional roughnesses. In this 
figure the abscissa is k/6*. It is seen that the three- 
dimensional trips are much more effective than the two- 
dimensional ones in promoting transition. That an in- 
crease in Mach number alleviates transition when in- 
duced by wires is also evident from the figure. 
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Summary 


The new flexural theory of elastic sandwich plates developed 
in two recent papers’ ? includes the effects of transverse shear 
deformations, rotatory and translatory inertias, and flexural and 
extensional rigidities, with no limitations imposed upon the mag- 
nitudes of the ratios between the thicknesses, material densities, 
and elastic constants of the core and faces. On the basis of the 
new sandwich plate theory and the exact elasticity theory, flexural 
vibrations of sandwich plates are investigated in this paper. 
Numerical results yielded by the two theories show excellent 
agreement with each other, and the new sandwich plate theory 
is seen to be good for a very wide frequency range that is of practi- 
cal interest. The importance of the various effects involved in - 
the theory is assessed. 


Symbols 
General: 
1 = (subscript) 1, 2, or 3, indicating the core, lower 


face, or upper face 
= rectangular coordinates in longitudinal and 
thickness directions, respectively 


Plate Quantities: 


hy = half-thickness of core 

hy = thickness of each face 

h = hi + he = total half-thickness of sandwich plate 
/ = length of plate in x direction 

Pi = density 

E; = Young’s modulus 

V; = Poisson’s ratio 

Nyy Mi = Lamé’s constants 

rh = 

Ty = 

= p2/pr 

= (1/m)[Ei/(1 — = — 

re = (1/m:)|E2/(1 — ve?)] = [2/11 — ry 


Vibration Quantities: 


w = circular frequency 
Q = (pi/p)w*h? 
= Vo/(1 + 
L = wavelength 
n = number of longitudinal half-waves in length of 
plate 
= 2zxh/L for infinite plates, or for finite 
plates 
é = d/h = wave number 
= 
4 In Sandwich Plate Theory: 
= bending moment 
Qc = transverse shearing force 
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Flexural Vibrations of Elastic Sandwich Plates' 
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Nz: = membrane force 

M, = Mn + + Mz 

Pz, Pz = functions of lateral loading in x and z directions, 
respectively 

u;, wu; = displacement components in x direction of zero 
and first order, respectively 

w; = displacement component in z direction of zero 
order 

w = transverse deflection of plate 

Yi, v2 = changes of slope of normal to middle plane of 


plate in core and faces, respectively 
W, vi, = amplitudes of w, 


«1 = E,/(1 — 

= E,/(1 — 

In = 2h,3/3 

= 2(he3/12) 

= 

= (2/3)(h3 — = Ino + + (h2/2)]*he 


kK, = ky = «x = shear coefficient (subscripts for identifying core 
and faces) 


In Elasticity Theory: 


Uj, Vir Wi = displacements 

e; = (Ou;/d0x) + (Ow;/dz) 

Cri = normal stresses 

Tyziy Tzxiy Tryi = Shearing stresses 

fis gi 

A Bi, given in (3.3) 

a;”, Bi? 

Cai?, 

ci? = (Ai + 2yi)/mi = 2[(1 — — 2%; )] 
a1, bi, dx, be \ defined in Eqs. (3.12) 

ay’, by’, as’, be’ 

A,, Ao, Az 

B,, Be, B;, Bs >» = defined in Eqs. (3.13) 

Ci, Cr f 

k 

X = ryull — (1/c2?)] { + rad) 


Introduction 


i Is a well-known fact that the effect of transverse 
shear deformation plays an important role in static 
problems of flexure of sandwich beams and plates.* * 
The shear effect is further known to be important in 
vibration studies of homogeneous beams and plates, 
as was pointed out by Timoshenko® and Mindlin.® 
The rotatory inertia was also considered by these 
authors, but its influence on the lowest branch of the fre- 
quency curve of an infinite plate was shown to be negli- 
gibly small.6 Since this is usually the only branch of in- 
terest to structural engineers, the impression has often 
been that the rotatory inertia is completely unimpor- 
tant, in spite of the fact that, if branches of the frequency 
curve higher than the lowest are to be brought into the 
theory of beams and plates, the rotatory inertia must be 
considered. As far as homogeneous beams and plates 
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are concerned, the higher branches cover rather high 
frequency ranges which have been so far of little in- 


terest to structural engineers. However, if and when 
information on such high frequencies is needed, it is then 
necessary to include these branches and thus the rota- 
tory inertia. Such is indeed the case with the vibra- 
tion of crystal plates, as has been investigated exten- 
sively and systematically by Mindlin.? 

In view of the above discussion on homogeneous 
beams and plates and the unknown situation with the 
vibration problem of sandwich beams and plates, a 
new sandwich plate theory has been developed,’ ? 
in which more factors are included than in previous 
theories. By assuming the transverse displacement 
to be uniform across the thickness of the plate and by 
assuming the displacements in the plane of the plate 
to be of linear variations over the thickness, with the 
slope in the faces taken to be different from that in the 
core, it is possible to include in the theory the transverse 
shear deformations in the core and faces, the rotatory 
and translatory inertias of the core and faces, the 
flexural rigidity of the core, and the flexural and ex- 
tensional rigidities of the faces. In addition, no limi- 
tations are imposed upon the magnitudes of the ratios 
between the thicknesses, material densities, and elastic 
constants of the core and faces of the sandwich. The 
theory is therefore actually applicable to symmetrically 
arranged three-layered plates in general. On the basis 
of the theory the free flexural vibrations of sandwich 
plates are investigated in the present paper, with par- 
ticular emphasis on ordinary sandwich plates that 
have thin but heavy and rigid faces. 

A brief summary of the new theory is first presented 
in Section (1). A frequency equation is next given in 
Section (2) for the infinite sandwich plate in plane 
strain. For ordinary sandwich plates the result may 
be simplified and is discussed in great detail. The im- 
portance of the various effects that are involved in the 
theory and the frequency equation is also assessed. 
To check the accuracy of the result obtained from the 
sandwich plate theory, a frequency equation of the in- 
finite plate is further derived in Section (3) from the 
exact elasticity theory. This result should also be of 
interest by itself. The sandwich plate theory, as any 
other approximate plate theory, has the advantage 
over the exact elasticity theory in that the former, 
although approximate, may be applied much more 
easily to finite plates. It is indicated in Section (4) 
how the result for infinite plates obtained from the 
plate theory may also be used immediately for plates 
that have finite lengths in one direction and have 
simply supported end sections. In Section (5) nu- 
merical examples are presented. Through one of these 
is shown the excellent agreement between the predic- 
tions of the sandwich plate and elasticity theories, 
thereby establishing the importance of the former 
theory. The numerical results also confirm the dis- 
cussion given in Section (2) concerning the importance 
of the various effects. Some important conclusions 
are finally drawn in Section (6). 
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(1) Résumé of Sandwich Plate Theory 


As in the previous papers,’ * the face layers of the 
sandwich plate are identical, and the middle plane of 
the plate is taken to be the xy plane. For the one- 
dimensional case of plane strain, only a unit length of 
the plate in the y direction need be considered. In the 
z direction, the thicknesses of the lower face, core, and 
upper face extend from —h to —h, —h, to ky, and hy 
to h, respectively. The thickness of the core is 
therefore 2/;, and that of each of the two faces is hy = 
h — h,. 

The stress-equations of motion of the theory are 


M,' — Q, + hpz = + 


po[(A® — hy) /3] + 
Ma’ + Nz’) On + hip, = (1,1) 


+ po [(h® — hy*)/2) X | 
+ te) + — tie) | 
+ pe = + pohe(ti2 + w;) 
in which primes indicate differentiation with respect to 
x and dots with respect to ¢. In terms of the plate dis- 
placements w, yi, and yY2, the displacement components 
are 


= YP 
= u;\)) = Yo 
= —u; = Yo) (1.2) 
w = w = w, = w 
Similarly, the plate stresses are 
Mn = yy’ 
=> M, ‘12) yp.’ + 
+ (1/2)hope" + (1/2)he] (1.2) 


On = + w’) 
= Qrs = + w’) 


Among these the bending moments and membrane 
forces have been expressed in a form somewhat different 
from that given in reference 1. In the present form 
it is easy to identify and as 
the flexural rigidity of the core and the flexural and 
extensional rigidities of each of the two faces, respec- 
tively, and yy, + (1/2)ho¥2 as the uniform extension 
of the faces. On substituting Eqs. (1.2) and (1.3) 
in Eqs. (1.1), the following displacement equations of 
motion are obtained: 


€1(241°/3) + 2eo(he®/12) po” + ) 
+ + — 

+ w’) — + w’) + | 

hp, = + + | 

+ + (1/2)hoPe) | 

€:(2h,3/3) + + (1, | ( 
+ w’) + = + 


(1.4) 


+(1 
+ w”) + + w”) + 
b: = 2(pili + 


the first two of which also appear in a different form 


) 
e | : 
7 
i 
1 
i 
6 
> 
> 


274 JOURNAL OF THE AERO/SPACE SCIENCE 


than those in reference 1. Similar to the rigidities, 
pi(2hy?/3), pe(h2*/12), pit, and poh. represent rotatory 
and translatory inertias of the core and faces. 

The equations of motion are subjected to the initial 
conditions of specifying the initial values of w, yi, Pe 
and their time derivatives, and to the boundary condi- 
tions of prescribing for the entire plate between + = 0 
/ one member of each of the two products 


Pr(ts)2—n- 


and prescribing for the end sections x = 
member of each of the three products 


[Mar + Nis — N22) 
+ — In( Nis — Ny) QO,w 


or, alternatively, the three products 


Min, [M2 + Ma — h(Na — N,2) |( 


and x = 
P:(W3)2=—n 


Qand x = /one 


v1), O,w 
(1.5) 


The method of determining the shear coefficient x 
for vibration problems was given in reference 2 and is 
analogous to Mindlin’s procedure for determining a 
similar coefficient for homogeneous plates. Briefly, 
the value of «x is determined by matching the lowest 
frequencies obtainable from the sandwich plate and 
elasticity theories for simple thickness-shear modes 
of the infinite sandwich plate. The frequency equa- 
tions for such modes according to the two theories were 
found to be, respectively, 


2 + + Q3 — 
Ki Ke pil pihy® pihy® 


+ 
= 


E 2 poh 
pily 


4+ (4%) 
pily Ko 


J (Fat) 
\ pilt,* 
4 
h,? pihy® 


polo 2 
pihy® kr \ way 


S—APRIL, 


1960 


(7 + — 12k[3r,r, + 14+ 


(rtn?/r,) + 36x? = 0 (1.6) 
tan f tan = 1/Vr,r, (1.7) 
For ordinary sandwich plates for which r,7,?/r, < 1, 


the equations may be simplified to yield approximately 
= [ror + (1/3) (1.8) 
ftanf = 1/rorn (1.9) 


The value of « is then determined by Eq. (1.8) after / 


is first solved from Eq. (1.9). It is easily seen that 


ll 


2/12 when 7,7 
1 r > 


the first of which is the case of a homogeneous plate. 


K 
xk = | approximately when 


(2) Vibrations of Infinite Plates Based on 
Sandwich Plate Theory 


The frequency equation may be derived by dropping 
the lateral-load terms involving p, and p, in Eqs. (1.4) 
and substituting in the results the following expressions 
of the plate displacements: 


w = W sin Ax/h | 
= V, cos Ax/h > (2.1) 
vo = W cos Ax/h 


where \ = 27h/L is a dimensionless quantity inversely 
proportional to the wavelength L. By following the 
usual procedure of setting equal to zero the deter- 
minant of the coefficients of W, Vi, and WY, and after 
some lengthy algebraic manipulation, we find 
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which is the frequency equation. The writing of the equation has been arranged in such a manner that all terms are 
dimensionless and that the effects of the various quantities may be readily detected. Among these quantities we 


have 


ply, = rotatory inertia of core 

polo = rotatory inertia of faces about their own middle planes 

polo, = rotatory inertia of faces about middle plane of sandwich plate (with faces taken to be concentrated at 
interfaces), actually contributed by translatory inertia po/2 of faces 

pele, = total rotatory inertia of faces, about own middle planes and about middle plane of sandwich plate, thus 
also involving effect of translatory inertia p2/» 

edi, = flexural rigidity of core 

éolo = flexural rigidity of faces about their own middle planes 


€:Jo, = flexural rigidity of faces about middle plane of sandwich plate (with faces taken to be concentrated at 
interfaces), actually contributed by extensional rigidity 2/2 of faces 
€:Jo2 = total flexural rigidity of faces, about own middle planes and about middle plane of sandwich plate, thus 


also involving effect of extensional rigidity eh 


The effect of any of these quantities may be suppressed by putting the quantity equal to zero in Eq. (2.2). The 
transverse shear effect in the core or faces is, however, suppressed by putting « or «2 equal to infinity in the equa- 
tion. In order to evaluate the relative magnitudes of the different terms, Eq. (2.2) is next written in terms of 
the various ratios and then divided throughout by 7,7,._ The result is 


2 8 1 
2(1 + ry)? (1 + 7,70) | — (2) - + + — - + 
2 8 1 ror,” 2 (2 
— (i + + i -- 
3 6 3 3 ) 


9 2 9 9 
(1 + 7)? | ve) + 27 (1 +% + + (2)(1 + (; n+ + 


8 1 ror,” 
— (2)(1 + r,rn) (; AP + 
Ke 6 
| + rn) +2 + 


l 2\4/8 1 rer,” 1 (2 


2 2 
fa + 7n)*(2)(2) rt (1 ++ At + 


8 2 
k (2) n+ ++) (20 += )| wh (23) 
3 6 Ke 3 4, 


Eqs. (2.2) and (2.3) have been arranged so that each of the terms may be traced back and forth immediately from 
one equation to the other. Being cubic in 2, they yield three branches of the frequency curve. 

Eqs. (2.2) and (2.3) are essentially the same as Eq. (47) of reference 1, in which the propagation of straight-crested 
waves in infinite sandwich plates was briefly discussed. The ratios 7,, 7,, and r, may have arbitrary values, and the 
equations are therefore applicable to any symmetrically arranged three-layered plate. For ordinary sandwich 
plates that have relatively thin but heavy and rigid faces, and for limited ranges of 2 and A’, these equations may be 
simplified. Thus, the value of 7,7, for such plates is usually of the order of unity and 7,/7, much smaller than unity. 
If values of Q are further limited to the order of unity, we have 


and allr,r,?/r, terms in Eq. (2.3) may be neglected. The equation then reduces to 
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+ ry)*(2)(2)(1 + + 


» 2 9 
(1 +7)? 13 + +r+ + = (2)(1 + (; n+ 


2 rer)? 
Ki Ke 3 6 


)+ 
Ky Ke 3 6 \ 


9 2 
(1 + r,)?(2)(2) + (1 + + 


8 1 l 5 ] 2 
3 /\6 3 ------- 


which is now of the second degree in 2, and from which 
only two branches of the frequency curve are obtain- 
able. 

Let us further assume that 


<1, <1 (2.6) 


the first being obviously true for most ordinary sand- 
wich plates. It follows from Eqs. (2.6) that, in the 
frequency equation, (a) all ror,*/r, terms may be 
neglected because the ratio r2/r, = 2/(1 — vw) is al- 
ways of the order of unity, (b) the term involving 7,7,°/6 
and i? is also negligibly small, and (c) so are the two 
r,°/3 terms. The frequency equation thus takes the 
simpler form 


(1/«)(1 + (1/3) + — 
+ + + + + 
rn) | + + rer |} + 
{[(1/3)r1 + rn) + [1/12«1(1 + X 
[(4/3)ri + rory = O (2.7) 


Since we usually have 


ret, > 1 


for ordinary sandwiches, the terms 1/3 + r,rn(1 + rp) 
and the 7; terms are actually negligible but have been 
kept for other future uses. 

Eq. (2.7) is important in that, although much simpli- 
fied, it will be seen to yield very accurate results for 
ordinary sandwich plates. The terms in the equation 
reflect the various effects that are still included in it. 
Tracing these terms back in Eq. (2.2), where they are 
underlined [as they are also in Eqs. (2.3) and (2.5)], 
we see that they involve the transverse shear effect in 
the core, the rotatory and translatory inertias of the 
core, the rotatory inertias of the faces about the middle 
plane of the sandwich plate (which is actually con- 
tributed by the translatory inertia of the faces), the 
flexural rigidity of the core (the effect of which is 
actually negligible), and the flexural and extensional 


rigidities of the faces. No longer included in Eq. 
(2.7) are the shear effect in the faces and the rotatory 
inertias of the faces about their own middle planes. 
The validity of Eq. (2.7) is restricted by the conditions 
set forth in Eqs. (2.4) and (2.6). However, as will be 
seen later, accurate results are obtainable from Eq. 
(2.7) even when some of these conditions are not quite 
fulfilled. 

A discussion of the effects of the various quantities 
that have been excluded from Eq. (2.7) as well as 
those that are included in the equation is in order. At 
the outset we mention that numerical calculations have 
indicated that the three roots of 2 of the frequency 
equation for a fixed value of \ are separated far apart 
in the case of ordinary sandwich plates. If we refer 
to these roots of 2 as the first, second, and third in the 
ascending order, which correspond to the first, second, 
and third branches of the frequency curve, the first 
root may be determined with very good accuracy by 
dividing the constant term by the coefficient of the Q 
term, the second by dividing the coefficient of the 2 
term by that of the 2 term, and the third by dividing 
the coefficient of the 2? term by that of the Q* term. 
It is thus seen that, by further referring to Eq. (2.2), 
the first and lowest root in the case of ordinary sandwich 
plates is strongly affected by the flexural and exten- 
sional rigidities, while the third and highest root is 
practically independent of these rigidities but depends 
much on the inertia quantities. Detailed discussion 
of the various effects now follows. 

The major influence of not including the shear effect 
in the faces or the rotatory inertias of the faces about 
their own middle planes is to reduce the degree of Q in 
the frequency equation from the third to second. The 
branch of the frequency curve thus eliminated, how- 
ever, is the third and corresponds to very high fre- 
quencies which lie beyond the practical range of applic- 
ability. It may therefore be disregarded without 
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affecting the usefulness of the remaining results. Fur- 
thermore, an additional branch of the frequency curve, 
whose cutoff frequency is that of the antisymmetric 
simple thickness-stretch mode, may exist between the 
second branch and the very high third branch given 
by Eq. (2.2) in the case of ordinary sandwiches. Since 
this additional branch is not contained in the present 
sandwich plate theory, the accuracy of the third 
branch may not be reliable, and its consideration is not 
warranted. In spite of the probable uselessness of the 
third branch for ordinary sandwich plates, however, it 
is not recommended that the shear effect in the faces 
or the rotatory inertias of the faces about their own 
middle planes be neglected so as to reduce the order 
of the equations of motion from the sixth to fourth in f¢. 
The reason is that this would not simplify the solution 
of any vibration problem of finite plates with edge 
conditions other than that of simply supported, as 
the order of the equations in x will be seen to be not 
reducible even for ordinary sandwiches unless we are 
interested only in rather low frequency ranges. On 
the other hand, to neglect the shear or rotatory inertia 
of the faces would destroy the generality of the equa- 
tions of the theory, which is actually applicable to any 
symmetrically arranged three-layered plate. 

If the shear effect in the core is further suppressed, 
Eq. (2.7) becomes 


+ + tn) + 
[(1/3) + + rm) — 
+ ren(1 + rm) = (2.8) 


The reason for having kept the terms 1/3 + 7,7, (1 + 
r,) in Eq. (2.7) is now clear. We see that the second 
branch of the frequency curve is no longer included in 
Eq. (2.8). Moreover, the accuracy of the equation for 
the first and lowest branch becomes questionable, since 
large ror, terms have been eliminated. In fact, nu- 
merical examples will show later that the results of 
Eq. (2.8) are simply not usable. 

We shall next investigate the influence of omitting 
from Eq. (2.7) the rotatory inertia of the core and 
that of the faces, both about the middle plane of the 
sandwich plate. As was indicated before, the latter 
is actually contributed by the translatory inertia of the 
faces. Eq. (2.7) in this case reduces to 


+ + + 
{((1/3)ra + + rm) + 
[1/12(1 + + X 
= (2.9) 


Although the second branch of the frequency curve is 
now also eliminated, as in the preceding case, the re- 
sults for the first and lowest branch given by Eq. (2.9) 
should be essentially the same as those given by Eq. 
(2.7), because the lowest frequency is given approxi- 
mately by the quotient between the constant term and 
the coefficient of the Q term, and because the terms 
1/3 + r,r(1 + ™) which no longer appear in Eq. 


(2.9) could have been omitted from Eq. (2.7) in the 
first place. 

Finally, if the flexural rigidities of the faces about 
their own middle planes are neglected, the \*® term drops 
out from Eq. (2.7), which then becomes 


(1/«i)(1 + + JQ? — 
+ + + + X 
(1 + + + + 
rotn |} + [(1/3)r1 + rorn(1 + = (2.10) 


Since the first two terms of Eq. (2.10) are identical with 
those of Eq. (2.7), the values of the higher root and thus 
the second branch of the frequency curve yielded by 
the two equations should be very much the same. On 
the other hand, because of the missing \* term in Eq. 
(2.10), the values of the lower root given by the two 
equations may differ considerably, especially for large 
values of \. The fairly large magnitude of the A* term 
is seen to be due to the joint effect of the flexural and 
extensional rigidities of the faces; the absence of either 
of the two eliminates the effects of both. While the 
parameter r2r,* (due to flexural rigidities) may be small 
by itself, the product between 72r,° and rer, (due to 
extensional rigidities) is usually a large number for 
ordinary sandwiches. Because of the relatively large 
influence on the lowest branch of the frequency curve, 
the term cannot be neglected, and the system of equa- 
tions of motion has to remain of the sixth order in x. 
This inherently complicates the vibration problem of 
sandwich plates. 

The complete frequency equation (2.2) or (2.3) may 
be reduced to Eq. (1.6) for simple thickness-shear 
modes by putting A = 0; namely, for infinitely long 
waves. Similarly, by taking A = 0 in Eq. (2.5) or 
2.7), Eq. (1.8) is obtained. In Eqs. (1.6) and (1.8) 
no distinction was made between x; and xo. 

The frequency equations (2.3), (2.5), and (2.7) for 
sandwich plates may further be reduced to the fre- 
quency equation for homogeneous plates composed of 


the core-layer material by putting 7, = 0. Thus, all 
three equations yield 
Q? [3x1 (Ca + = 0 (2.11) 


A similar frequency equation for a homogeneous plate 
composed of the face material may also be obtained by 
setting 7, = © in Eq. (2.3) but not in Eq. (2.5) nor 
Eq. (2.7), the result being 
[(r,/7,)2]? — {3x2 + + (r2/r,) X 
[(r,/r,)2] + ko(T2, r,)r4 = 0 (2.12) 
Since the one-dimensional case of the plate is being 
considered here, Eqs. (2.11) and (2.12) correspond 
essentially to the frequency equation derivable from 
the Timoshenko beam equation.’ If, on the other 
hand, we set 7, = Oorr, = © in Eq. (2.8), there results 
22 ‘=0 
or (3 + A*)r,Q — rdAt = O 
which corresponds to the frequency equation of a 
classical homogeneous beam. 


| 
| 
| 
1S 
| 
| 


(3) Vibrations of Infinite Plates Based on 
Elasticity Theory 


The frequency equation for the flexural vibration of 
an infinite sandwich plate in plane strain will now be 
derived from the exact elasticity theory. The dis- 
placement equations of motion are 


+ (Ai + (Oe;/Ox) = 

+ (Ai + mi) (Oe,/02) = 
where ¢e; = (Ou;/Ox) + (Ow,/Oz) and 7 = 1, 2, or 3 for 
the three layers of the sandwich as before. The general 
solution of Eqs. (3.1) is the same as that for a single- 
layered plate, which is given, for instance, in reference 
7. The final form of the solution is as follows: 


= (Ef; gi’) cos x et 

(fi’ + sin 

= + — + 2g," | x 
siné xe 

= — BP)f, + 2Eg,’| sin x 

Oy, = —Ai(a? + sin — x onl 

+ (#2 — B2)gi] cos x 


Since the case of plane strain is considered, the dis- 
placement v, and the shear stresses 7,,;, 7,2; are zero. 
In Eqs. (3.2) the prime indicates differentiation with 
respect to z, — is related to \ by A = &h, and f; and g; 
are functions of z of the form 


(3.1) 


Font 


Teri 


bo 


fi A; sin a + B; COS 

gi = C; sin 8,2 + D; cos 
with af +2 =w/ca?, BE +E = 
= (Ai + 2ui)/ Pir 


(3.3) 
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A;, B;, C;, and D; in f; and g; are arbitrary constants. 
The boundary conditions at the free surfaces of the 
sandwich plate are 


=0 —f) (3.4) 


Continuation of the displacements and stresses at the 
interfaces between the layers requires that 


= (2 = —h,) (3.5) 


Us, = W3, On = O23, Tel = Tae (2 = 


Since the flexural motion is antisymmetric with respect 
to z, f; and g; are taken in the following form: 


fi = A; sin | 
= Ao sin acs + Bo cos avs | 

3.6 
= dD, cos Biz ( 


=F Cy sin Boz + cos Boz 
3 


With /; and g; so chosen, only the six conditions given 
in Eqs. (3.4) and (3.5) need be considered. We sub- 
stitute Eqs. (3.6) in Eqs. (3.2) and the results in turn in 
Eqs. (3.4) and (3.5), to yield six homogeneous equa- 
tions, linear in the six constants .1;, D;, Ao, Bo, Co, and 
D.. The determinant of the coefficients of the six 
constants is then set equal to zero in the usual manner. 
After some lengthy algebraic manipulation the following 
frequency equation is obtained : 


— 1)(@* — (2g? — — cos ayy sin + 
— Be”) — — Bi*)] — Bi?) — 7, — Bo”) ] sin ayhy cos + 
{ — r,(&? — cos sin + — Bi?) — — sin cos Bin} X 
{ 4£% COS cos Boh. — — sin sin Bohs} _ 
{4(r, — cos sin + [27,é? — B,*) | sin cos Byhy} x 
4£a082 sin sin Beh, — — Bo”)? cos ashy cos Bohs} + 


+ B,7)(E? + Bo?) cos ay cos X 


| cos sin Boh, + — Bo”)? sin cos Boho} — 


roi + + Be) sin ayy sin X 


sin ashy cos + (£2 — Be”)? cos sin Boho} = (3:7) 


Reduction of Eq. (3.7) to the frequency equations 
for some simpler special cases may readily be carried 
out. By taking h; = 0 or h, = r, = 0, Eq. (3.7) be- 
comes 


a8; cos a,h; sin + 
(#2 — 6,7)? sin cosB A; = 0 (7 = 1, 2) 


which is the well-known Rayleigh-Lamb frequency 
equation for flexural vibration of a homogeneous plate 
and has been discussed in detail by Mindlin.”.* For 
very long waves, we take é = 0, and Eq. (3.7) may be 
shown to yield 


Bi = he = O 
,( a2) tan + tan = (3.8) 
r,(B2/8;) tan tan Beh. — 1 = 0 


Since we now have 


w/Ceis Bi = 


the first two of Eqs. (3.8) simply yield w = 0. The 
last two of Eqs. (3.8) are the frequency equations for 
the simple thickness-stretch and simple thickness-shear 
modes of sandwich plates, respectively, the latter of 
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which has been discussed in reference 2 and is the same 
as Eq. (1.7), but with different notations. 
For further discussion, Eq. (3.7) will now be rear- 


ranged to appear in a somewhat simpler form. For this 

purpose we write, by virtue of Eqs. (3.3), 

= — & 
= [1/1 +n) IV @/a2) (3.9) 
= + (P/a?) — 1 = ada 

where 

c® = + = 2101 — — 

r = (3.10) 


= 
| 


+ — (P/a2) 


ay 


I being thus proportional to the square of the phase 


velocity. It follows from Eq. (3.9) that 


a; = hy )ay = 1é(1 + = 


a’ = (1 = v1 (r C1") 


where 


Ay = a’b,’, 
B, = Ar, — 1) 
B; = r,(1 + — 2 

Ci = — — be’) 
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In a similar manner we rewrite 8,/, azhs, Bohs, Bi, ae, 


and @.. Tosummarize, we have 
ah, = tray a, = ta,’ 
(3.11) 
avhy = a = 1£d2 | 
Bohs = 1X bo Bo = 
ay’ 14(1 + rp) 
a, = = h 
+ ln \! a" 
(3.12 
| rr 1+r, 
1+ Pr, C2” Th 
r ] + ln 
= — W1- by’ = by - 
l + Tn \ J 


By making use of the relations in Eqs. (3.11) and the 
following additional notations: 


1+ 6," 

B, = (1 + b,’*) 

By = r,(1 + be”) — (1 + by") 
Cy = — — be"*) 


(3.13) 


the frequency equation (3.7) may be rewritten in the following form: 

tanh Aq, [4.124;B.By — A(4B,? + cosh Aas cosh Abo + 
+ A37B,?) sinh Ade sinh Ay tanh Nb, [4.40.1;B,B; Ao(4B;3? + A;37B,?) x 
cosh cosh Abe + (442°B,? + sinh sinh Ade] + C,[4A2 cosh sinh Abo — 


A;? sinh cosh + C. tanh Aq; tanh Ab, [4-42 sinh Aa» cosh — 13? cosh day sinh Ady] = 0 


It should be mentioned that the replacement of the 
trigonometric functions by the hyperbolic functions in 
the frequency equation does not serve any special 
purpose. 

Eq. (3.14) may best be used for determining \ for 
given I. The values of c,” are first computed from the 
Poisson ratios of the face and core materials according 
to Eqs. (3.10). With given values of 7,, r,, 7, of the 
plate and a chosen value of I’, the values of a, bi, ay’, 
b;’,... may now be computed from Eqs. (3.12), and 
A,, By, GQ, ... in turn from Egs. (3.13). A is then the 
only unknown in Eq. (3.14), which may be solved by 
trial and error. A further step will yield 2, once \ 
becomes known. 

It has been found difficult to apply Eq. (3.14) directly 
to the calculation of the lowest branch of the frequency 
curve for ordinary sandwich plates, where the ranges of 
\ and T are limited, since it is not practically possible 
to obtain accurate results. For such cases, the equa- 
tion is rearranged by introducing series expansions of 
the hyperbolic functions such as 


cosh Adz = 1 + (A*a2?/2) + (A4a24/24) + 
+... 
sinh Adz = Ade[1 + + 
+... ] 


(3.15) 


(3.14) 


Inside each square bracket in Eq. (3.14) terms are then 
arranged in ascending powers of \*. The coefficients 
of most of the terms are complicated functions of de, de, 
ly, A3, By, Bo, Bz, and By; and, in their final form, the 
coefficients may be arranged to appear as polynomials 
of powers of k = (I/2)(r,/r,). The constant terms 
and k& terms in the coefficients all vanish; the poly- 
nomials start with k* terms. Since 7, is usually much 
larger than r, for ordinary sandwich plates, k can be 
much smaller than unity if we limit ourselves to values 
of T small enough, and the polynomials of k may be 
truncated. Take for example the coefficient of \* term 
in the second bracket in Eq. (3.14), which is 
(44 °B,? + A;°B3*)debe 
16.42[%/(1 + rm) — 1) [1 — — 
{2r,2[1 — (1/c2?)] + 2r,(1/e2?) — 
[1 + (1/c2?) + 7,2h4) = 
32 Ask*[rn/(1 + ra) — — 1) — 
Most of the other coefficients are much more involved 
than this one shown here, but they may be similarer 
truncated. In this manner, Eq. (3.14) becomes, with 
some additional small terms involving k furthyl 
neglected, 
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xe{e hd 4 = ( 2 a) (1 


where 
X = — (1/2?) + m) IA}? (8.17) 
Eq. (3.16) was thus deduced from Eq. (3.14) under 
the assumptions 
k = (P/2)(r,/r,) <1 
{ + r,)]\}? = small and <1 


\ (3.18) 


: which, for a given ordinary sandwich plate, limit the 
ai use of the equation to small values of T and A; namely, 
to the lower portion of the lowest branch of the fre- 
quency curve. For Eq. (3.16) to be applicable to large 
values of I’, the polynomials of powers of k, which con- 
stitute the coefficients in the \? series, should not be 
truncated. Similarly, to extend the use of Eq. (3.16) 
for large \, more terms should be kept in the series ex- 
pansions in Eqs. (3.15), which will result in more terms 
in the two series in the X? term of Eq. (3.16). How- 
ever, both of these processes are extremely tedious. 
For small enough \ values, on the other hand, the 
Ks third and even the second terms in each of the two 
i series in Eq. (3.16) may not be needed. The smallness 
of the quantity { [7,/(1 + rn) |r} 2 determines the num- 
ber of terms to be kept in the series, which is the im- 
plication of the second assumption given in Eqs. (3.18). 

The procedure of computation by the use of Eq. 
(3.16) begins with assuming a pair of values of T and 
\; for example, those obtained from the sandwich plate 
theory may be used. With assumed values of I and 
A, Eq. (3.16) may be solved for X. \ may in turn be 
calculated from Eq. (3.17), and the result checked 
against its assumed value. When the calculated and 
assumed values of \ are not close enough, this procedure 
should be repeated; for instance, by taking the calcu- 


4 lated \ value to be the new assumed value. The nu- 
: merical example to be given later tells us that a second 
: calculation is usually not necessary when assumed 

values of T and X are taken from the results of the 


i sandwich plate theory. 


(4) Vibrations of Finite Plates Based on 
Sandwich Plate Theory 


For a plate that has a finite length / in the x direction 
and is simply supported at the end sections x = 0) and 

= /, the boundary conditions may be taken to be, 
according to (1.5), 


it 
a‘b’ tanh 1 + + — = (; + 108 (; +, 
6X {tanh (+1) - 1] — a,;’b,;’ tanh -1] - 


a’ 
a) +h 


2 
3 jtanh day E 4 — a,’b;’ tanh \d; — a 


1 
(ann Ad, tank AD 
[rn/(1 + - 


Tr 
= 0 (8.16) 
4\l+n, 


~ 
° 6 
& 
| 
E | 
€0.(2.10) 
c | 
8 | 
A 
02 o6 10 12 


w =0 


M, = €1(2h,3/3) + po’ 
Qeoho(hy + he /2) (hiy’ + /2) = 0 
+ — Nis — = + 
+ + howe’ /2) = 0 


These conditions are satisfied at x = 0, / if the plate 
displacements w, ¥1, and y2 given in Eqs. (2.1) are again 
used, with the exception that \ is now taken to be 
nth/l, n being the number of half-waves contained in 
the length / of the plate. The frequency equations de- 
rived in Section (2) are therefore immediately appli- 
cable to such simply supported finite plates with this 
new interpretation of 2. 


(5) Numerical Examples 


In the following examples the core and face materials 
of the sandwich plate are assumed to be cellular cellu- 
lose acetate and aluminum respectively, for which r, = 
34.4, r, = 1683, 7; = 2.20, and re = 4790. Two values 
of r, are chosen—namely, 1/2 and 1/10. Since r,r, > 2, 
the shear coefficient « is taken equal to 1 according to 
Eq. (1.10). The calculations were limited to the ranges 
VQ < 2.0 and < 1.4, and two branches of the fre- 
quency curve were included. 

Results of the calculations for the case r, = 1/2 
have been plotted in Fig. 1. Since Eqs. (2.4) were 


satisfied, calculations based on the sandwich plate 
theory were made according to Eq. (2.5) instead of 
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(2.3). In addition, Eq. (2.7) was made use of, even 
though Eqs. (2.6) were not quite satisfied for the right- 
hand portion of the lower branch. As expected, re- 
sults for the lower branch from the two equations (2.5) 
and (2.7) deviated somewhat from each other as \ 
increased. On the other hand, since \ was small for the 
entire range covered by the upper branch in Fig. 1, 
results for this branch given by the two equations were 
so close that the difference could not be made visible 
in the figure. Eq. (2.8), which does not include any 
shear effect, and Eq. (2.9), which does not include the 
rotatory inertia of the core and faces about the middle 
plane of the sandwich plate, were next employed in 
the calculation, and both yielded only the lower branch 
of the frequency curve. As had been expected, the 
result of Eq. (2.8) deviated greatly from that of Eq. 
(2.7), but the frequency values given by Eq. (2.9) were 
only very slightly lower than those given by Eq. (2.7). 
The difference between the results of Eqs. (2.9) and 
(2.7) become greater as \ increased, but within the range 
of calculation the maximum difference at \ = 1.4 was 
only about 1 per cent and could hardly be indicated in 
the figure. Eq. (2.10), which does not include the 
flexural rigidities of the faces about their own middle 
planes, was further used to compute the frequency 
curve. The result for the higher branch was essentially 
identical with that of Eq. (2.7), but the lower branch 
obtained was considerably lower than that given by 
the latter equation, the deviation becoming much 
more severe as \ increased. 

The frequency equations (3.14) and (3.16) derived 
from the elasticity theory were finally used for calcu- 
lating the two (lowest) branches of the frequency curve 
for the case 7, = 1/2, the former for the higher branch 
and the latter for the lower one. When Eq. (3.16) was 
used, assumed values of T and A were taken from the 
results of Eq. (2.5) of the sandwich plate theory. The 
value of k decreases rather rapidly with \. For \ equal 
to 1.4, 1.0, and 0.8, the corresponding values of k 
are approximately (0.026, 0.013, and 0.009. The condi- 
tion set forth by the first of Eqs. (3.18) was therefore 
fulfilled rather satisfactorily for the range covered by 
the lower branch in Fig. 1, and the use of Eq. (3.16) 
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02 04 06 os L2 


Fic. 2. = 1/10. 


20 


w 


Fic. 3. rm, = 0. 


was justified. Since the highest value of [r,/(1 + 
r,) JA in this range was 1.400/3 or 0.467, the third term 
in each of the two series in Eq. (3.16) was small enough 
to be neglected. Results for the lower branch ob- 
tained from Eq. (3.16) were almost the same as those 
from Eq. (2.7), and a separate curve corresponding to 
the former was not inserted in Fig. 1. Such a separate 
curve would be only very slightly higher than that 
given by Eq. (2.7), with deviation between the two 
curves beginning to be noticeable at about A = 1.0. 
Similarly, results for the higher branch obtained from 
Eq. (3.14) could be considered for all practical purposes 
to be identical with those of the sandwich plate theory. 

Results for the case r, = 1/10 were shown in Fig. 2. 
Since the requirements stated in Eqs. (2.4) and (2.6) 
were essentially satisfied, Eq. (2.7) was used. The two 
branches of the frequency curve calculated are ex- 
pected to be in close agreement with the results of the 
exact elasticity theory just as in the previous case. 
Also the same as in the previous case are the results 
obtained for the present case from Eqs. (2.8) and 
(2.9), which both yielded only the lower branch— 
namely, while the result of Eq. (2.8) differs consider- 
ably from that of Eq. (2.7), the result of Eq. (2.9) is 
essentially the same as the latter. Similarly, for the 
present case, Eq. (2.10) again yielded a lower first branch 
than Eq. (2.7). Although the deviation between the 
two is smaller for this than for the previous case, it is 
nevertheless of a substantial amount except when d is 
smaller than, say, 0.6. The upper branch given by 
Eq. (2.10) is again practically identical with that 
given by Eq. (2.7). 

In addition to the cases r, = 1/2 and 1/10 for the 
sandwich plate, calculations were made for the case 
r, = 0 for a homogeneous plate composed of the core 
material. Eq. (2.11), which includes the effects of 
transverse shear and rotatory inertia, gave two 
branches, while Eq. (2.13), which does not include these 
effects, yielded only the lower branch. The results 
have been plotted in Fig. 3. The two lower branches 
are seen to deviate from each other, just as in the case 
of a sandwich plate, but the deviation is much smaller. 
The second branch is seen to start off at a much higher 
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cutoff frequency than the sandwich cases, and it is also 
much flatter. 

A word should be added here if the frequency curves 
are to be used for simply supported finite plates. In- 
stead of a continuous variation of \ for the infinite 
plate, \ takes now only discrete values. The curves 
in Fig. 1, and similarly in the other two figures, are 
meaningful only at such discrete values of \ when 
applied to finite plates. The points A, B, C, D,..., 


Ai, Bi, .. . in Fig. 1 represent some such values of \ 
for the case of h/] = 1/20. For h// = 1/10, only 
every other point such as A, C,..., Aj, .. . should be 
used. 


(6) Discussion and Conclusion 


Some important and interesting conclusions may 
now be drawn from the results of the above examples. 
The results in the case 7, = 1/2 for both branches of 
the frequency curve plotted in Fig. 1 show excellent 
agreement between the predictions of the elasticity 
and sandwich plate theories. There is slight deviation 
between the theories for the right-hand portion of the 
lower branch, and the results of Eq. (3.16) (elasticity 
theory) are found to be closer to those of Eq. (2.7) than 
to those of Eq. (2.5) (both of the sandwich plate 
theory). One source of error, as mentioned before, was 
that the first assumption in Eq. (3.18) involved in de- 
riving Eq. (3.16) was only approximately fulfilled for 
this portion of the lower branch. The accuracy of 
Eq. (3.16) may always be improved, but additional 
calculation would be rather lengthy and _ tedious. 
Since the difference among the various curves is not 
large at all, this is considered unnecessary. We thus 
conclude that Eq. (2.7) has included all the essential 
elements for ordinary sandwich plates and is satisfac- 
tory even for values of r, as large as 1/2. As 7, de- 
creases, the agreement between the results of the sand- 
wich plate theory and the exact elasticity theory is 
expected to become closer. 

The numerical results have further confirmed the 
discussion in Section (2) on the importance of the vari- 
ous effects that are included in Eq. (2.7) but are sup- 
pressed one by one, respectively, in Eqs. (2.8) to (2.10). 
Thus, the very large deviation between the results of 
Eqs. (2.7) and (2.8) has established the importance of 
the shear effect in the core (just as in static flexural 
problems of sandwich plates). The shear effect is in 
fact far more important for sandwich than for homo- 
geneous plates in vibration studies. 

Next, we observe that the cut-off frequency for the 
second branch—namely, the frequency of the simple 
thickness-shear mode—is much lower and closer to the 
first branch in the sandwich than in the homogeneous 
case. (If the sandwich plate in the above examples 
has a total thickness of 3/10 in., the cut-off frequency 
is 6,700 cps for r, = 1/2 and 10,600 eps for 7, = 1/10, 
both of which are within the audible range.) This 
frequency and, together with it, the second branch of 
the frequency curve are thus much more important in 


the sandwich case. Consider, for instance, for 7, = 
1/2 a simply supported finite sandwich plate with a 
ratio h/] = 1/10. The gravest frequency is naturally 
given by point C on the lower branch in Fig. 1. The 
second gravest frequency, however, is not given by the 
next higher point E on the same branch but is equal 
to the cutoff frequency of the upper branch indicated 
by point A;, which has a lower value than that given 
by E. Since it would not be possible to include the 
upper branch of the frequency curve in the plate 
theory without considering the rotatory inertias of the 
core and faces about the middle plane of the sand- 
wich plate (the latter being actually contributed by the 
translatory inertia of the faces), their importance is 
obvious. As was pointed out before, the negligible 
effect of the rotatory inertia on the lowest branch of 
the frequency curve should not be interpreted indis- 
criminately as an indication of its being unimportant. 

Finally, a factor of unexpected importance to ordi- 
nary sandwich plates, even for values of 7, as small as 
1/10, has been found to consist in the joint effect of the 
flexural rigidities of the faces about their own middle 
planes and their extensional rigidities, which vanishes 
if either of the two types of rigidities is disregarded. 
For a wavelength equal to 7 times the total depth of 
the plate (A = 1) in the case 7, = 1/10, the frequency 
is 12.5 per cent too low if this joint effect is not taken 
into account. For larger values of 7, and for shorter 
wavelength (larger \) the inaccuracy increases very 
rapidly. Except for plates with very thin faces and for 
sufficiently low frequency ranges, therefore, the flexural 
rigidities of the faces about their own middle planes 
must be included even though they may be small by 
themselves. As a consequence, the system of equations 
of motion will have to remain of the sixth order in x 
and cannot be reduced, which at once complicates the 
vibration problem of finite sandwich plates with bound- 
ary conditions other than that of simply supported. 
For plates with very thin faces and for sufficiently low 
frequency ranges, some simplification of the vibration 
analysis of sandwich plates appears possible but will be 
treated separately. On the other hand, the present 
analysis is seen to be good both for rather large values of 
r, and for a very wide frequency range that is of practical 
interest. The range may cover up to, say, 20,000 cps 
or higher, as long as it is not too much higher than the 
first thickness-shear frequency. 
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Magnetohydrodynamic Effects on Exothermal 
Waves' 


R. A. GROSS,* W. CHINITZ,** anv T. J. RIVLIN*** 
Fairchild Engine and Airplane Corporation 


Summary 


Elementary magnetohydrodynamic effects on exothermal, 
plane, supersonic waves are examined. One-dimensional steady 
flow with heat addition in the presence of a plane magnetic field 
transverse to the direction of motion is studied. Basic equations 
are derived and solved numerically for several ideal conditions 
which correspond to gaseous detonations and thermonuclear 
shocks. The magnetic field acts on the flow in a way which 
tends to counteract the heat-addition effects. The singular 
solution still exhibits the classical Chapman-Jouguet conditions. 
Some properties of high-temperature gas mixtures are presented 
including the electrical conductivity of air shocks and fuel-air 
detonations. Experimental implications of the analysis and 
numerical results are discussed. 


Symbols 
A(T) = electron - atom, molecule, orion collision cross-section 
B = magnetic induction vector 
7 = propagation velocity of a small disturbance 
Diy. = diffusion coefficient 
E = internal energy 
E = electric-field vector 
H = magnetic-intensity vector 
A = H?*/p 
J = current-density vector 
k = Boltzmann constant 
L = length 
m = particle mass 
mM, = electron mass 
M = Mach number 
n = number density of electrons 
p = pressure 
Q = heat added 
Q = Q/RT 
ri2 = distance between particle centers at collision 
R = gas constant 
Rp = ratio of magnetic force to dynamic force 
Ry = ratio of induced magnetic field to applied magnetic 
field 
Ss = entropy 
t = time 
T = absolute temperature 
u = ratio of final to initial velocity 
V = velocity vector 
a = mole fraction of ionized constituents 
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= ratio of specific heats 
= electron charge 
magnetic permeability 
density 

= electrical conductivity 
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ll 


(1) Introduction 


era is a subject in which one 
studies the motion of electrically conductive 
fluids subject to magnetic fields. This area of study, 
once the sole domain of the astrophysicist, has recently 
attracted the interest of investigators in many diverse 
fields. 

We seek here to investigate the question of what role 
magnetohydrodynamics may play in exothermal flow, 
such as gaseous detonations or thermonuclear shocks. 
The analysis considers only macroscopic phenomena 
with emphasis upon the behavior of the fluids. Ideal- 
ized cases are studied to formulate the important phys- 
ical parameters and to gain an understanding of the 
effects that might be expected. In particular, one- 
dimensional steady flow with heat addition is analyzed 
for the case in which the magnetic field is plane and 
transverse to the direction of motion. The equations 
that describe such a flow are solved numerically and 
graphs are presented showing the kinds of effects to be 
anticipated for stoichiometric detonations and fusion- 
type shock waves. It is shown that in the presence of 
a magnetic field the analogue of the Chapman-Jouguet 
detonation (singular-point solution) has a product 
Mach number still equal to one and the product en- 
tropy is stationary. 

Calculations were performed to obtain the electrical 
conductivity of the product flow of fuel-air strong deto- 
nations and air shocks, and these results are shown up 
to a Mach number of twenty (about 12,000°R.). 
Experimental implications are discussed. The effect 
of the Earth’s magnetic field on previous combustion 
studies is shown to be negligible. 


(2) Fundamental Equations 
(@) Frozen Flux Relation 


The behavior of a conducting fluid is thought to be 
governed by the equations of fluid dynamics, Maxwell’s 
equations, and Ohm’s law. The latter equations, in 
the rationalized mks system are 


(OB/dt) + curl E = 0 (la) 
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H 


Fic. 1. Schematic illustration of the flow. 


curl (B/u) = J (1b) 
div B = 0 (Ic) 
where B = wH 
and Ohm’s law is 
J =o(E+V XB) (2) 


where we have assumed that the displacement currents 
are negligible. Let us restrict our attention to the case 
where the plasma motion is steady and one-dimen- 
sional and uw is constant. Thus, in particular, 


V = [V(x), 0, 0] (3a) 
where 
V = V; = constant —-ocx<x 
V = V(x) aS 
V = V2 = constant 
and B = (0, 0, R(x) ] (3b) 
where 
B = B, = constant —-xax<x<x 
B = B(x) Sx 
B = B, = constant x < 


Thus, in the steady regions upstream and downstream 
of the wave all quantities are considered constant and 
designated with a subscript of 1 or 2. Therefore, in 
regions ] and 2, 


VXE=40, 


Therefore, from Eq. (1b) the current J must be zero 
in these regions of no change. Ohm’s law, Eq. (2), 
therefore implies that 


Ey = —V, x B, = ViB, 
E, = x B, = 


VxB=0 


But, since the flow field is steady for all values of x, 
Eq. (1a) gives 


VXE=0 


and therefore the electric field E has only a y component 


which is everywhere equal to a constant. Hence, 
Fe 
or VB, = (4) 


Eq. (4) is known as the frozen flux relationship. It is 


1960 


important to note that this relationship, for the steady 
one-dimensional case treated here, was arrived at 
independent of the value of the electrical conductivity. 
This fact was kindly brought to our attention by Dr. 
M. Mitchner. The frozen flux relationship is usually 
derived by assuming infinite electrical conductivity 
and employing Eqs. (la) and (2), but here, for steady 
one-dimensional flow which contains regions of constant 
properties, the frozen flux relation is valid for gases 
with finite conductivity. 

In the following parts of this paper we will employ 
nonrationalized, Gaussian notation so as to be in accord 
with previous related papers, such as references 2 and 
14. 


(b) Propagation of Small Disturbances 


Consider now the question of the propagation rate of 
a small plane disturbance through a conducting, com- 
pressible, homogeneous medium subject to a magnetic 
field. We assume that the flow and magnetic field are 
given by Eqs. (3a) and (3b). Conservation of mass is 
expressed by 


d(pV) = 0 
or Vdp + pdV = 0 (5) 


where p is the density and V is the velocity. Conserva- 
tion of momentum is expressed by 


dp + pVdV + (HdH/4n) = 0 (6) 


where / is the pressure and H7 the magnetic field 
strength. Combining Eqs. (5) and (6) we obtain the 
propagation velocity of a small disturbance, c, given by 


c? = (dp/dp) + (H/4m)(dH/dp) (7) 
Employing Eqs. (4), (5), and (7) gives 
c? = (dp/dp) + (H?/4mp) (7a) 


The first term of Eq. (7a) is the square of the speed of 
sound in any elastic fluid, while the second term is the 
contribution of the magnetic field and is identical to 
the square of the Alfvén-wave speed. 

Conservation of energy in this one-dimensional, 
steady, adiabatic case is expressed by 


VdV + [yR/(y — 1)]dT + 
(HdH/2ep) — = 0 (8) 


where 7 is the specific-heat ratio, 7 is the absolute tem- 
perature, and R the gas constant. Employing Eq. 
(7a), Eq. (8), and the ideal-gas equation of state gives 


c? = yRT + (H?/4rp) (9) 


Thus, in one dimension, the square of the speed of 
propagation of small disturbances in a conducting fluid 
is the sum of the square of the usual gasdynamic 
speed of sound and the square of the Alfvén-wave 
propagation rate. In several dimensions the problem 
is more complex and results in several propagation 
speeds. 
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(c) Conservation Equations 


The equations which express conservation of mass, 
momentum, and energy for a compressible, electrically 
conducting fluid have been derived in detail by Cowl- 
ing,| deHoffmann and Teller,? Burgers,* and Pai.* 
For cases in which the frozen-flux relation applies, the 
magnetic field acts like a hydrostatic pressure /7?/Sx 
together with a tension //*/4m along the magnetic lines 
of force. 

For one-dimensional steady flow, such as indicated 
in Fig. 1, Eqs. (5), (6), and (8) can be integrated when 
the electrical resistance is negligible. Thus, conserva- 
tion of mass is expréssed by 


= (10) 


Conservation of momentum is expressed by 
pi + + = po + p2V2 + (11) 


Conservation of energy with heat addition is expressed 
by 


(Vi?/2) + [yRTi/(y — 1)] + + = 


(V22/2) + [yRT:/(y — 1)] + (12) 


where Q is the thermal energy added to the flow by com- 
bustion, fission, fusion, or other forms of energy release. 
Eqs. (4a) and (10) combined give a familiar form of the 
frozen-flux relationship, 


= pe (13) 


(3) Solution of the Equations 
(a) Solution of the Conservation Equations 


To determine the effect of the magnetic field on the 
flow parameters, we combine Eqs. (10) through (13), 
together with the equation of state, to obtain a cubic 
equation for the velocity ratio (final to initial velocity) 
in terms of the initial conditions only. The result is 


(Vo/V;)? — xy(Vo/V;)? + + w =0 (14) 
where 
2 
+ 
y = 14+ + + 
s = [(y — 1)/2]M2[1 + + 
((y — + 
— 1)/v(Q/RT)) + 1 
and 
M? = 


yRT + 


Eq. (14) was solved approximately by the Newton- 
Raphson'* method for various values of the dimen- 
sionless parameters 


y, ff = Q = Q/RT; 


Graphs of the solutions V2/ V, and 72/7; plotted against 
M, are shown in Figs. 2 and 3. 

It is particularly interesting to examine the singular 
solution of Eq. (14)—the analog of the Chapman- 
Jouguet solution in the classical case. The coefficients 
x, ¥, 2, w, are all positive for any realistically given 
parameters and initial conditions. Therefore, accord- 
ing to Descartes’s Rule of Signs, the number of positive 
roots of Eq. (14) is either two or zero. Hence, in the 
case in which there are positive roots, Eq. (14) has 
two positive roots and one negative root. The singular 
case is the one in which the two positive roots coalesce. 

We compute the singular values by determining the 
minimal value of J/, that leads to a solution of Eq. (14). 
This we do by solving Eq. (14) for M, in terms of u = 
V2/ Vi, and minimizing the resulting function of u in the 
range 
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The singular solutions are indicated by circles in Figs. 
2and 3. 

Having found V2/V; numerically, other properties 
may be readily computed. For example, it readily 
follows that 

p2/pr = Vi/V2 (15) 

= 1 + (1 — (Vi/V2)?] + 
+ — (V2/Vi)] (16) 
T2/T, = (p2/pi)(V2/ Vi) (17) 

In the classical case we have M, = | at the Chapman- 
Jouguet point. Computation of an example (Q = 
100, y = 1.667, 7 = 10) leads to My = 1.0002 at the 
singular point. Indeed, as we show now, the presence 
of magnetic fields does not alter the conclusion that at 


the singular point ds = 0 and Mz = 1 
Eq. (12) may be written as 


+ + (Vi?/2) + = 
Ex + + (V2?/2) + (p2/p2) (18) 


where /; and /» are the internal energies at states 1 
and 2, respectively, and £;, includes the heat of reaction 
Q. Putting 


E* = E + (H?/8xp) (19) 
and p* = p + (H?/8r) (20) 


and substituting in Eq. (18) we obtain the conservation 
of energy equation in the form 


E,* + + = 
+ (V2?/2) + (p2*/p2) (21) 


Eq. (11) may now be written as 
pi* + = po* + poV2? (22) 


Upon eliminating the velocities among Eqs. (21), 
(22), and (10) we obtain the Hugoniot relationship 


E,* — E\* + (1/2)(p2* + pi*) X 
[(1/p2) (1/1) = 0 (23) 


Now considering state 2 to be variable and differenti- 
ating Eq. (23), we obtain 


dE* + (1/2)(p* + pi*)d(1/p) + 
(1/2)[(1/p) — (1/m) = (24) 


It is a consequence of Eqs. (10) and (22) that 


(p* — pi*)/[(1/p) — = 
—(pV)? = —(piVi)? = constant (25) 


In the (p*, 1/p) plane, segments [i.e., lines satisfying 
Eq. (25)] from (p;*, 1/p;) that meet Eq. (23) (for vari- 
able state 2) do so in either 2 or 1 points (corresponding 
to solutions shown in Figs. 2 and 3). When the seg- 
ment is tangent to the curve we have the singular condi- 
tion, the analog of the C-J condition. When this 
happens we have at the point of tangency 


dp* = [(p* — pi*)/(1/p) — (1/p-)|d(1/p) (26) 
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The second law of thermodynamics can be written as 
Tds* = dE* + p*d(1/p) (27) 
and Eqs. (24) and (26) when combined with Eq. (27) 
yield 
ds* = 0 


at the point in question. In view of Eqs. (19), (20), 
and (13), the first line in Eq. (27) implies 


Tds* = Tds = dE + pd(1/p) (28) 
or ds* = ds (29) 


At the singular point (p.*, 1/p2) we have on the one 
hand 


dp*/d(1/p) = (30) 
in view of Eqs. (10), (25), and (26), and on the other 
hand 


dp*/d(1/p) = [dp/d(1/p)|s + 


(H2/4m) [dH /d(1/p)] 


— p2*(dp/dp)s — 
(dH /dp) 


where differentiation along the Hugoniot curve at the 
singular point can be replaced by differentiation along 
the adiabatic curve at the singular point by Eq. (29). 
Egs. (30) and (31) together give 


= (dp/dp)s + (H2/4n)(dH/dp) = 


according to Eq. (7). Hence we arrive at the im- 
portant conclusion that at the singular point, 


= V2? =] 


(b) Discussion of Results 


From the solution of Eq. (14) shown in Figs. 2 and 
3, we can observe the following general properties: 

(1) For fixed Q, an increase in /7 drives the singular 
solution toward AJ, = 1. 

(2) For fixed /7, an increase in Q drives the singular 
solution away from M, = 1. 

(3) For fixea M, and Q, an increase in 7 drives the 
two solutions (if they exist) toward their limits of 1 and 

(4) For fixed M, and H, an increase in Q, drives the 
two solutions (if they exist) away from their limits of 1 
and [(y — 1)/(y + 1)]. 

(5) The singular solution approaches 1 as M, ap- 
proaches 1. 

The magnetic field effects upon the flow tend to op- 
pose the influence of the heat addition effects on the 
flow. The general character of these solutions exhibit 
the same characteristics as simple plane detonations; 
at a given Mach number there is either no solution, 
one solution (singular C-J solution), or two solutions 
(strong and weak detonations). The transition from 
the detonation case with no magnetic field to one with 
a magnetic field takes place continuously, so that a 
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small magnetic field acts like a perturbation to the 
nonmagnetic combustion problem. 

The values chosen for Figs. 2 and 3 correspond to 
physically interesting cases. A typical stoichiometric 
hydrocarbon-oxygen mixture reacts with 0 ~ 100. 
The interesting case of thermonuclear fusion represents 
Q = 1,000. A magnetic field strength of 5,000 gauss 
is equivalent to a pressure of approximately 1 atmos- 
phere. Field strengths of 10,000 gauss are readily 
achievable in the laboratory today. Similarly, pres- 
sures of the order of 10~* atm. can be easily achieved 
using present vacuum techniques. Consequently [7 
can be varied from 0 to 10°. The value of y = 5/3, 
corresponds to a fully ionized gas in which the particle 
velocities are random and only three degrees of freedom 
are present. 


(4) Some Properties of High Temperature Gas 
Mixtures 


Because of the one-dimensional and steady-flow 
assumptions we demonstrated in Section (2) that the 
frozen flux effects are independent of the value 
of the electrical conductivity, so long as it is non- 
zero. There are many conditions where one cannot 
approach one-dimensional steady flow and then elec- 
tron current flow results and the electrical conductivity 
assumes great importance. We were interested to 
learn whether combustion may contribute to signifi- 
cantly increasing the electrical conductivity of fluids 
in chemical equilibrium under conditions that are now 
of engineering importance. To assess the properties 
of high-temperature gas mixtures we have carried out 
some computations on a digital computer. These 
numerical results permit an evaluation of the electrical 
conductivity, chemical composition, pressure, temper- 
ature, etc., of gas mixtures that have been brought 
to high temperature by chemical and aerodynamic 
means. The general reaction considered is 


C,H, + + + — + 
+ + + BsCO + + 
B;N + BsO + + + 
Bull, + BeA + + + 
BysNO2 + Bie( NOt + e-) 


where x and y are integers which determine the fuel, 
a; are the moles of reactants and 6; are moles of prod- 
ucts. The products are in thermodynamic equilibrium 
and the problem is solved for steady one-dimensional 
flow. Details of the techniques employed in these 
computations are contained in reference 5. 

The product temperatures of air shocks and some 
stoichiometric fuel-air strong detonations are given 
in Fig. 4. The reactant initial conditions were fixed 
at 537°R. and 10-* atm. It is readily seen that disso- 
ciation is significant above 2,000°R. and that above 
Mach 8.5 (5,500°R.) the addition of a chemical fuel 
does little to increase the product temperature. In 
fact, it may reduce the product temperature, the energy 
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Fic. 4. Product gas temperature vs. shock or detonation Mach 
number. 


going into the dissociation of molecules such as //,0, 
CO», etc. 

In an effort to gain an appreciation of the effects 
of ionization on the electrical conductivity of the 
product gases, we have utilized the previous numerical 
results to determine ¢. The equation used to deter- 
mine the conductivity was that proposed by Chapman 
and Cowling® for slightly ionized gases composed of 
rigid elastic spherical molecules: 


= 
(ne®/RT) \(3/Snry2) [RT + 


where 
n = number density electrons 
e = electronic charge 
k = Boltzman constant 
7 = absolute temperature 
Dy = diffusion coefficient 
‘2 = distance between particle centers at collision 
m = mass of particle 


If we introduce the collision cross section A(7) defined 


A(T) = 


and assume that the mass of the electron is negligible 
compared to that of the molecule, we obtain finally 


o = [1/A(T)]-10° 
or 


= 5.16 X mbhos/em. 


where 
a = mole fraction of ions 
m, = mass of electron (gms) 
A(T) = electron - atom collision cross section (cm.?) 
T = absolute temperature (°R.) 
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The value of the collision cross section A(7) was 
calculated by adding the cross sections (weighted with 
the number density) of each molecule and atom which 
contributed significantly to the total cross section. 
Values for this parameter were obtained from references 
7 and 8 and ranged from 9 K 10~-" cm.? to 4 K 107" 
em.?. 

In Fig. 5 is shown the electrical conductivity of the 
product gases for air shocks and strong detonations up 
to Mach 20. These theoretical calculations are based 
upon the assumption of thermodynamic equilibrium 
and the only ionized molecule considered is NO*. 
Nitric oxide has a relatively low ionization potential 
(9.5 volts) and is the major source of ions at the tem- 
peratures considered. For example, for air at 7) = 
12,415°R. and p = 0.586 atm. 


Mole Fraction 
NO+ 
Nt 
O+ 
4% 10-* 
O.+ 


To produce large magnetic effects on the fluid it is 
desirable that the electrical conductivity be large. The 
electrical conductivity is governed primarily by the 
mole fraction of ionized constituents and hence the 
temperature and chemical composition of the gas 
mixture. The alkali elements have the lowest ioni- 
zation potentials (Cs = 3.89 volts, Na = 5.138 volts, 
K = 4.34 volts) and hence might be used to ‘‘seed”’ 
the mixture to increase the concentration of electrons. 
Fig. 5 shows the effect of seeding air with cesium. The 
concentrations shown are one molecule of cesium to 10! 
and 10° molecules of the mixture. The Saha equation 
was used to determine the equilibrium constants for 
the alkali metals. 

In addition to the equilibrium ions considered pre- 
viously, detonations create appreciable quantities of 
ions resulting from the chemical reaction. These 
ions, resulting from the so called chemi-ionization 
process, have been discussed by Calcote.’ In subsonic 
flames, measurements have been made which indicate a 
conductivity of ¢ ~ 3 X 10-* mhos/cm. in the com- 
bustion zone. This is far in excess of what one would 
calculate using thermal-equilibrium considerations. 
Thus, the chemi-ionization process in a detonation 
wave may be of assistance in producing magnetohy- 
drodynamic effects. 


(5) Experimental Implications 


Several dimensionless parameters help to give some 
perspective as to when large magnetohydrodynamic 
effects might be expected in real gas mixtures. Resler 
and Sears” have discussed some of these parameters. 
The ratio of the electric body force to the dynamic 
pressure, designated by Rg, is 


Rp = off? V1/(1/2) pV? 
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This is physically similar to a Reynolds number, and 
when Rz > 1 we expect significant magnetohydrody- 
namic effects. The value of R, is shown in Fig. 6 for 
some detonation and air shock conditions. For Mach 
numbers in excess of about 15, for the particular condi- 
tions chosen (JJ = 104 gauss, Z = 5 in., convenient 
laboratory values), significant magnetohydrodynamic 
effects may be expected. 

The ratio of the induced magnetic field to the applied 
magnetic field is 


which has been called by some persons a magnetic 
Reynolds number, although physically Rg, is more 
analogous. Values of Ry are shown in Fig. 7 and for 
the conditions considered it can be seen that the in- 
duced magnetic field is small compared to the applied 
magnetic field. 

The validity of assuming that the properties of a 
fluid may be specified by its macroscopic nature, with- 
out delving into its kinetic nature, was raised by 
Bullard.'' He states that the macroscopic theory is 
valid if the radius of curvature of the electron orbit 
paths (cyclotron or gyro radius) is large compared to 
the mean free path. That is 


H <mVnA(T)/e 
From kinetic theory, 
H 
or H < 3.78 X 10°[pA(T)/VT] 


where p is in atmospheres, 7 is in °R., A(7) is in em.’, 
and #/ is in gauss. The regions where the gyro radius 
exceeds the mean free path are shown in Fig. 8. 

However, when ¢@ is large, Spitzer'® states that the 
macroscopic motion of an electrically conducting fluid 
is not sensitive to the ratio of collision frequency to 
cyclotron frequency. This view is taken by Chapman 
and Cowling,® who state that even if the interval be- 
tween two consecutive encounters is very long, a uni- 
form steady state (independent of kinetic effects) is 
possible. Consequently, the experimentalist must 
examine carefully his particular objective and decide 
whether the ratio of gyro radius to mean free path is 
important. 

It is clear that in all combustion research to date the 
effect of the Earth’s magnetic field has been negligible 
because of the small value of 7, the Earth’s field being 
approximately 0.2 gauss. This results from the fact 
that the previous solutions exhibit the property that 
small values of H, act only as a perturbation to the 
problem where H = 0. The fine work on one-dimen- 
sional shocks in shock tubes by Mitchner' and Patrick 
and Brogan” indicate that one-dimensional, adiabatic, 
steady flow of a compressible fluid in a magnetic field 
can be adequately predicted. Consequently, it seems 
hopeful that our similar approach to the analysis of 
the nonadiabatic case with its corresponding different 
class of solutions may be verified by experiment. The 
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unbounded one-dimensional case given here is experi- 
mentally unrealizable but may be approximated by a 
central core of fluid buried in a large flow field. Such 
approximations as this may some day come into per- 
spective in much the same way as did inviscid fluid 
theory when it was finally coupled with boundary-layer 
theory. 
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Shear Lag Solutions for Sheet-Stringer Panels 
by Means of a Hydrodynamic Analogy 


MARTIN GOLAND* 
Southwest Research Institute 


Summary 


An analogy is established between the stress flow in flat sheet- 
stringer panels and the plane potential flow of an incompressible 
fluid. The sheet-stringer panels are presumed to be “‘shear lag”’ 
structures carrying essentially longitudinal loads (parallel to the 
stringers), and the effects of lateral stresses and deflections are 
ignored. 

It is shown that the hydrodynamic analogy consists of a direct 
relationship between longitudinal stress and the longitudinally 
directed component of fluid velocity. Shear stress bears a direct 
relationship to the laterally directed component of fluid velocity. 
In transferring from the sheet-stringer panel to the equivalent 
flow, an affine transformation of coordinates, and, hence, of the 
panel boundaries, must be made. 

A discussion is given of the corresponding boundary conditions 
for the stress flow and the fluid flow. A boundary free of normal 
stress is equivalent to flow past a solid wall. It is also shown how 
rigidly constrained boundaries can be dealt with in relatively 
simple fashion. Problems of reinforcement around cutouts in 
the sheet-stringer panel are also discussed. 

The use of the analogy is demonstrated in several examples, 
including the stress concentration around elliptic cutouts, with 
free and constrained boundaries, in large panels under uniform 
tension. Also studied is the case of a concentrated longitudinal 
force applied at the center of a large panel. 


(1) Introduction 


7 PROBLEMS dealt with here relate to the stress 
analysis of flat plate structures composed of a 
cover sheet stiffened longitudinally by stringers and 
laterally by ribs. The applied loads are presumed to 
be primarily longitudinal and in the plane of the plate. 

For structures of this kind, it is often possible to 
approximately predict the main features of the stress 
flow by means of a “shear lag’”’ analysis. The essence 
of such an analysis lies in neglecting the influence of 
lateral stresses and deflections on the longitudinal and 
shear stress flow through the structure. The results of 
a shear lag analysis, when properly applied, give reli- 
able estimates of the principal stressing conditions— 
i.e., the maximum longitudinal stress in the stringers 
and their associated effective width of cover sheet and 
the maximum shear stress in the cover sheet. 

In this paper, a relatively simple analogy is estab- 
lished between the flow of longitudinal and shear 
stresses in a plate-like ‘‘shear lag” structure, and the 
potential flow of an incompressible fluid. The analogy 
is useful for obtaining solutions applicable to various 
shear lag problems, particularly those involving cutouts 
in the plate, and also affords a clear physical visuali- 
zation of the stress flow. 


Received February 2, 1959. 
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(2) Establishing the Analogy 


A schematic sketch of a typical shear lag structure is 
shown in Fig. l(a). The longitudinal loads are applied 
in the x-direction, at the boundaries of the plate; and 
the structure shown has a rectangular cutout which 
interrupts two stringers and removes six shear panels. 
The rib spacing is L, and the stringers are spaced a dis- 
tance ) apart. As shown in Fig. 1(b), the cover sheet 
thickness is ¢, and the cross-sectional area of a stringer 
is A,. For the moment, it will be assumed that all 
stringers and shear panels have identical properties 
over the entire plate. 

Considering now a typical bay of the structure, it is 
seen that the effective area for the transmission of longi- 
tudinal stress is the stringer area 4, plus the effective 
width of sheet on each side of the stringer. If the cover 
sheet is unbuckled, the effective area per bay is A = 
A, + bt. It will be convenient to define an “‘effective’’ 
thickness for longitudinal stress transmission as 


t, = A/b (1) 


Shear stresses are transmitted entirely through the 
cover sheet, and it is assumed that the sheet remains 
fully effective in shear. 

Fig. 2 shows a free-body diagram for a typical bay 
consisting of a length of stringer extending between 
two ribs plus the sheet extending (1/2)) on each side. 
Let o be the mean of the longitudinal stresses acting 
over the width of the bay and let r be the mean of the 
shear stress acting over the length of the bay. Equilib- 
rium of the bay is then expressed by the difference 
equation 
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+ + + + 
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(tr + Ay) Ax 


(tt) Ax 


Fic. 2. Equilibrium of a typical bay in the plate structure. 


[A(at,)/Ax] + [A(rt)/Ay] = 0 (2) 


The equations of elasticity for the plate element can 
be written as the difference equations 


Au/Ax = 0/E, Au/Ay = (3) 


where w is the longitudinal displacement, / is the modu- 
lus of elasticity in the x-direction, and G is the shear 
modulus of the sheet material. 

While the discussion can be continued on the basis 
of the difference equations, Eqs. (2) and (3), it is more 
convenient to generalize the problem concept as fol- 
lows: imagine a continuous plate having the fictitious 
characteristic of a thickness /, in resisting longitudinal 
stress o, and a thickness ¢ in resisting shear stress r. 
The lateral stresses and deflections of the plate are 
ignored. The equations of equilibrium and elasticity. 
comparable to Eqs. (2) and (3), are then 


(Qot,/dx) + (drt/dy) = 0 (4) 
Qu/dx = o/E, du/dy = 7/G (5) 


Eq. (5) at once leads to the compatibility condition 
governing the strains in the plate; using the concepts 
of generalized plane stress this can be written as 


0°u/Oxdy = (1/Et,)(Oot,/Oy) = (1/Gt)(Ort/Ox) (6) 


Now define the dimensionless distance variables 


x’ =x/L, y’ = y/b (7) 
and also define 

S = at,/L, q = rt/b (8) 

Then, Eqs. (4) and (6) become 
(0S/Ox’) + (0q/dy’) = 0 (9) 
0S/dy’ = B(dg/dx’) (10) 

where B is the shear lag factor 
B = (E/G)(t,/t)(b?/L*) = (E/G)(A/bt)(b°/L*) (11) 


Finally, make one further coordinate transformation as 
x'/VB = x/LVB, =y/b (12) 


and, on the basis of Eq. (10), define the potential func- 
tion $(&, 7) as 


S = VB(d¢/dt), = (13) 
Eq. (9) then becomes 
= (0%/dE) + (0%/dn?) = 0 (14) 


1960 


The shear lag problem has thus been reduced to the 
solution of Laplace’s equation in the coordinates (é, 7). 
These coordinates are related to the plate coordinates 
(x, y) by the affine transformation given in Eq. (12). 

Moreover, if ¢ is visualized as a velocity potential in 
the plane problem of the hydrodynamic flow of an in- 
compressible fluid, it is then seen that the horizontal 
flow velocity at a given point can be interpreted as 
the quantity S/ VB, while the vertical velocity is 
simply equal to the shear flow gq. 

Thus, in the interior of the plate, the analogy is estab- 
lished. 


(3) Discussion of the Boundary Conditions 


In the actual structure, two conditions involving 
either the applied stresses or deflections along the plate 
boundaries must be specified in order to make the prob- 
lem unique. Thus, at a free edge, both the normal and 
tangential stresses must be zero; for a completely con- 
strained edge, the deflections are zero, and so forth. 

The approximate analogy under discussion here does 
not permit the satisfaction of two conditions at the plate 
boundaries; rather, the unique specification of the har- 
monic solution suggested by Eq. (14) allows only the 
specification of either ¢ or its normal derivative along 
the various portions of the plate boundaries in the 
(~, ») plane. Hence, some physical difficulty should 
be anticipated when arriving at a suitable approxima- 
tion for the analogy boundary condition. 

If the deflections are specified at the plate bound- 
aries, the difficulties evaporate because of the neglect 
of the influence of lateral deflections. Thus, com- 
bining Eqs. (5) and (13), it is readily shown that 


u = (b*/Gt)d (15) 


so that lines of constant potential ¢ represent lines of 
constant longitudinal deflection. In particular, along 
a fully constrained boundary ¢ = u = 0. 

Consider now boundaries where the state of stress 
is specified. The simplest manner of arriving at the 
appropriate boundary condition is through the vari- 
ational derivation of the equations of equilibrium, as 
outlined by Love.'! It must be kept in mind that the 
present analysis is based on satisfaction of the equation 
of equilibrium in the longitudinal (x-) direction only. 
The application of Eq. (5), p. 167 in Love, then leads 
to our Eq. (2). The corresponding surface condition 
given by Eq. (6) of Love is then 


(OW/0e,,) cos (x, v) + (OW/de,,) cos (y, v) = X, (16) 


where W can be interpreted as the strain energy function 
per unit area of plate; ¢,, and e,, are the longitudinal 
and shear strains, respectively; X, is the surface load 
per unit length of boundary acting across the tangent 
plane v in the x-direction; and cos (x, v) and cos (y, ») 
are the direction cosines of the outward normal to the 
plane v. 
For the present analysis, it is readily shown that 
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OW = Et,(0u/dx) = LW B (d¢/dE) (17) 
OW/de,, = Gt(Qu/dy) = b(d¢/dn) (18) 


Next consider Fig. 3{a), which shows an element of 
the plate boundary in the (x, y) plane. It is clear 
from the sketch that cos (x, v) = cos @, cos (y, v) = 
Fig. 3(b) shows the corresponding element of 
Then 


sin 6. 
the transformed boundary in the (£, 7) plane. 
it is readily seen that 


tan y = —(dt/dn) = 7 
— (dx/dy)(dt/dx)(dy/dn) = (b/LVB) tan@ (19) 


and 


cos y = dn/dt = (dy/ds)(dn/dy)(ds/dt) = . 
cos sin? @ + (20) 


Using Eqs. (16)—(20), the boundary condition for 
stresses then becomes 


(0¢/0E)cos y + (0¢/On) sin y = 
(X,/LV B)[(1/B)(b/L)? sin? + cos? (21) 


It is to be noted that the left side of Eq. (21) is merely 
the fluid velocity normal to the boundary in the (£, 7) 
plane. 

From Eq. (21), it is at once clear that a stress-free 
boundary—i.e., a boundary free of normal stress, in the 
(x, y) plane corresponds to flow past a solid wall in the 
(é, n) plane. 

Problems of considerable practical importance relate 
to the treatment of structures with a reinforced area. 
Thus, in the main region of the plate the shear lag 
factor may have the value B,, while in the reinforced 
area the value of the shear lag factor may be B:. In 
principle, the hydrodynamic analogy is capable of ap- 
plication to such problems. Two separate flow planes 
must now be considered for the two regions of plate. 
Along the common boundary (in the (x, y) plane) be- 
tween the two regions, the conditions to be satisfied 
are: 

(a) The deflections u in both regions must be the 
same. The corresponding condition imposed on the 
velocity potentials for the two regions at the common 
boundary is given by Eq. (15). 

(b) The surface forces X, representing the action 
of each region on the other must be equal in magnitude 
and opposite in sign. The fluid velocities normal to 
the common boundary in each flow plane must then 
be related through Eq. (21). 

For the general reinforcement problem, it is found 
that analytic difficulties preclude the practical use of 
the hydrodynamic analogy. In the case when B, = 
B:, however, the analogy can be conveniently applied. 
(Note that the shear lag factor is simply the ratio of 
the effectiveness of the plate in resisting longitudinal 
vs. shear stresses. If, for example, both ¢, and ¢ are 
increased over a reinforced area but their ratio is kept 
the same as in the unreinforced plate region, the value 
of B will remain constant throughout the structure.) 

An illustration of the use of the analogy for a rein- 


forcement of this type around an elliptic cutout is given 
later. 


(4) Some Results for Elliptic Cutouts 


As a first illustration of the use of the analogy, con- 
sider an elliptic cutout in a large stiffened panel under 
uniform tension. The uniform tension acts in the x- 
direction and is of amount 7 per unit length of y- 
coordinate. The axes of the elliptic cutout lie along 
the (x, y) axes and are of length (2x, 2yo), respectively. 

Consider first the case when the elliptic cutout is 
unreinforced. The elliptic boundary in the (x, y) 
plane then transforms into an ellipse in the (£, 7) plane 
with the axes 


= = 2y0/b (22) 


It is now necessary to introduce some results for the 
plane hydrodynamic potential flow in the (Y, Y) plane 
past an elliptic boundary. For this purpose, introduce 


the curvilinear elliptic coordinates (a, 8), such that 
X =ccoshacos®, Y=csinhasin8 (23) 


The contours a = constant are confocal ellipses with 
their major axes along the X-axis. For the ellipse de- 
fined by a = apo, let 2a) be the major axis and 2b) the 


Xd Normal 


\ 


Tangent Plane uv 


dy 


Fic. 3(a). Plate boundary in the (x, y) plane. 
"Normal 
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Fic. 3(b). Equivalent boundary in the (£, 7) plane. 
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minor axis. Then, ad) = cosh ap, = c sinh ae, and 
= Ao? by?. 

Suppose that the flow in regions remote from the 
elliptic boundary a» consists of a uniform flow with ve- 
locity V in the Y-direction—i.e., normal to the major 
axis of the ellipse. Then the velocity potential and 
stream function for the flow past the boundary ap are 
given by* 

(24) 


(25) 


g= V} Y + aodre~ sin B} 
y= Vv} —X + * cos B} 
where 
= [(ao + b0)/(ao — bo) 


If the flow is normal to the minor axis of the boundary 
ajp—i.e., in the X-direction, then 


V{X + bodre-* cos B} 
V{Y — bde~* sin 


(26) 


(27) 


Now return to a consideration of the analogy be- 
tween the shear lag structure and the fluid flow. Sup- 
pose, first, that the elliptic cutout has a stress-free 
boundary. The hydrodynamic model is then the flow 
of a uniform stream in the £-direction past a solid 
elliptic boundary with axes 2, 2n0 as previously de- 
scribed. In the event 7 > &, the model must be 
interpreted with the use of Eq. (24). For & > m, 
Eq. (26) is appropriate. It can then be shown that the 
maximum concentration of longitudinal stress occurs 
at the points & = 0, » = +m. The maximum stress 
concentration factor is given by 


K, = 1 + (m0/é0) (28) 
and, using Eqs. (11) and (22), 
Ky = 1 + (y0/x0) [(E/G)(A/bt) (29) 


Eq. (29) can be applied to the case of an unstiffened 
plate of uniform thickness by setting 4 = bt and 
E/G = 2.6 (for a Poisson’s ratio of 0.3). While the 
shear lag model is not strictly applicable to the case 
of an unstiffened plate, it is interesting to compare the 
predictions from Eq. (29) with those from the exact 
theory of elasticity. From Eq. (29), 


K, = ] + 1.6(y0/xo) 


Results from exact theory are identical in form, but 
the factor 1.6 is replaced by 2.0. 

Consider next the circumstance when the boundary 
of the cutout is rigidly constrained. Then u = ¢ = 0 
everywhere along the boundary. But y = 0 along 
the ellipse ap for both Eqs. (25) and (27). It is at once 
clear that these equations are applicable to the present 
problem when they are interpreted as_ velocity 
potentials. 


(30) 


* In hydrodynamics, it is conventional to define the (x, y) com- 
ponents of the fluid velocity in terms of a velocity potential ¢, 
such that u = —¢gz,7 = — gy. 
negative signs are omitted. 
made for the stream function. 


In the present discussion, the 
A similar change in treatment is 
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Carrying the argument further along these lines, it 
can be shown that for the rigidly constrained boundary 
the maximum longitudinal stress concentration occurs 
at the points & = +, 7 = 0. The stress concentra- 
tion factor is of amount 


K, = 1 + (&/no) (31) 
i.e. 


K, = 1 + (x0/¥0) [(G/E) (bt/A)]'” (32) 


For a circular cutout in a homogeneous plate («») = 
yo, A = bt, E/G = 2.6), Eq. (32) gives K, = 1.62. 
Results from the exact theory of elasticity give the 
value K, = 1.52. 

Finally, consider the circumstance when the elliptic 
cutout has a stress-free boundary, and the plate is rein- 
forced over a region surrounding the cutout. Let the 
subscript 1 refer to the main region of plate outside of 
the reinforced area, and let subscript 2 denote the rein- 
forced region. Then, the reinforcement will be as- 
sumed to be of the special type such that B,; = Bz, 
I, = Ls, and 6; = bo. This type of reinforcement can 
be achieved by increasing both ¢, and ¢ in proportion 
over the reinforced area; suppose 


= to/t = 6 


Two flow regions must now be dealt with in the 
hydrodynamic analogy, corresponding to the regions 
1 and 2. However, since B; = Bs, etc., the common 
boundary between the two regions appears as a geo- 
metrically identical curve in both the (&, m) and 
(£2, n2) planes. 

For mathematical convenience, it will be assumed 
that the common boundary between the two regions, 
as it appears in the (£, m2) plane, is a confocal ellipse 
to the ellipse representing the cutout boundary trans- 
formed to the (&, m2) plane. To illustrate the method 
of solution, consider the case when these ellipses have 
their major axes along the £-axis. Let &, & be the 
major semi-axes for cutout and common boundary, 
respectively, and let 4, 7, be the corresponding minor 
semi-axes. 

Now introduce curvilinear elliptic coordinates (a, 8) 
such that 


£=ccoshacos8, »=csinhasin8 (33) 


where c? = &? — mo”. Then region 2 is defined by 
ao < a < a,, where tanh ap = m/f and tanh a, = 


ne/é-. Region 1 is defined by a > a,. 
For the velocity potentials in the two regions, take 
¢: = V[c cosh a + vmde~*] cos B (34) 
= V2[e cosh a + mAe~*] cos (35) 


where 


= + — m0) 


and where the values of v and V2 are to be determined 
through the conditions at the common boundary. 
It is to be noted that in regions remote from the 


re 


Ww 


wl 


whe 


Nov 


ct 
ol 
di 

| 

| 
fre 

y= 

tio 
| 
me 
cul 

: avi 
rei 
ter 

(5) 
2 ] 
lag 
lati 
pla 
loa 
4 cen 
of 
rest 
q 

cast 
I; 
int 
alor 
| long 
app 
of s 
in tl 
F 
flow 


) 


SHEAR LAG SOLUTIONS FOR 


cutout (a —~ ©), the flow in the (£), ;) plane consists 
of a uniform velocity (tension) of amount V in the é- 
direction. Since 0¢g:/O0a@ = 0 at the boundary of the 
cutout (a = a»), this boundary is stress-free. 

Now, at the common boundary (a = a,), Eq. (15) 
requires that 


= = (36) 
while Eq. (21) requires that 
= 0¢;/0a (37) 


These two conditions enable the evaluation of v and V> 
from the equations 


+ nom) = V6(E + vqom) (38) 
Vo(ne — nom) = — (39) 
where 
(& — — m0) 


With the values of v and V2 known, the stress distribu- 
tion throughout the structure is at once determinate. 

Since the present discussion is concerned mainly with 
methodology, and, because of the somewhat imprac- 
tical profile chosen for the common boundary, the cal- 
culations are not carried further. Analytic tools are 
available, however, to permit the study of practical 
reinforcement profiles following the same general pat- 
tern for the solution. 


(5) Concentrated Longitudinal Force Applied 
at the Center of a Large Plate 


In a relatively recent and excellent paper on shear 
lag problems, McComb? deduces a method for calcu- 
lating the stresses around rectangular cutouts in flat 
plate structures. As part of his method, he uses cer- 
tain influence coefficients representing the effects of unit 
loadings. One such unit Joading consists of a con- 
centrated longitudinal force applied at the intersection 
of a rib and a stringer. It is interesting to note the 
results afforded by the hydrodynamic analogy for this 
case. 

In the hydrodynamic analogy, a concentrated force 
applied to the center of a large plate can be interpreted 
in terms of a source outflow from an infinite wall placed 
along the y-axis. To see this, note first that if the 
longitudinal force is directed along the x-axis and is 
applied at the origin (x = 0, y = 0); then, by virtue 
of symmetry, the y-axis is free of normal stress. Hence, 
in the flow plane, the y-axis becomes a solid wall. 

For flow in the half-plane — > 0, place a source out- 
flow at the origin (£ = 0,7 = 0). The velocity poten- 
tial in the presence of the solid wall is then 


Inr 


(40) 


where 


+ 


r 


Now, using Egs. (20) and (21) and integrating around 
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a semicircle at the origin with the vanishing radius e, 


f (O¢/On) «dé = We On)dt = 
LV B)(cos y/cos (41) 


where 0¢/0n is the fluid velocity normal to the contour. 
But 
(cos y/cos 6)dt = (dn/dy)ds = (1/b)ds 


and so Eq. (41) can be written as 
(1/2)0 = (1/oLV B) X,ds = Pw/bLWVB (42) 


where Po is the total force applied at the origin to the 
half-plane x > 0. For a unit fotal force applied to the 


entire plate, then Po = 1/2. Hence, 
= (LV B)— 


In his analysis, McComb uses the scheme shown in 
Fig. l(a) for designating the location of stringer-rib 
intersections and shear panels by means of the sub- 
scripts 7, 7. Consider now the total longitudinal force 
associated with intersection 7j, comprising the sum 
of the forces in the stringer and the effective width of 
cover sheet on each side. In terms of the present 
solution this can be expressed as (coordinates x’ and 
y’ from Eq. (7) are used) 


j+1/2 
Pi; b( ot.) x" =i dy’ = 
j 


—1/2 


j+1/2 
LbV B dn (43) 
j 


But 
bLV ve = (WB/2m)[i/(i2 + Bn?)] (44) 
and so 


Py = ( x 
[i+ (1/2)] — tan-(WB/i)[j — (1/2)]} (45) 


Eq. (45) replaces the very much more complex integral 
expression given by McComb in his Eq. (A18). 

A comparison of the numerical values calculated by 
McComb for various P;; with those resulting from 
Eq. (45) shows differences on the order of 10 per cent 
in regions close to the applied load. In view of the 
different approximations involved in the two solutions 
closer agreement is not to be expected. 


Concluding Remarks 


The main emphasis in this paper is on the demonstra- 
tion of the analogy between stress flow and fluid flow. 
The application of the analogy to the treatment of a 
variety of problems of practical importance must come 
in later studies. 

Perhaps the principal significance of relating the 
stress and fluid problems is to show clearly how the 

(Continued on page 303) 
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A Class of Optimum Trajectory Problems 
in Gravitational Fields 


E. W. GRAHAM 


Douglas Aircraft Company, Inc. 


(1) Summary 


Problems involving the transfer of a rocket vehicle from one 
point to another with minimum fuel expenditure are considered. 
The rocket vehicle is assumed to operate in a gravitational field 
without atmospheric resistance. The transfer time and the 
terminal velocity vectors are specified. 

It is already known that, if the magnitude and direction of the 
gravitational force are everywhere constant (a uniform field), 
then the fuel required is a minimum if impulses are applied only 
at the terminals. This result is here extended to a class of prob- 
lems in centrally directed fields. A dimensionless parameter P 
involving the strength of the gravitational field and the pre- 
scribed transfer time is used to measure the importance of the 
nonuniformity of the gravitational field. When P = 0 the field 
is effectively uniform and the optimum procedure is to use termi- 
nal impulses with a coasting trajectory between the terminals. 
When P is greater than zero but less than 1.0 it is still true, for a 
restricted class of terminal velocity vectors, that the use of 
terminal impulses with a-coasting trajectory between terminals is 
optimum. 

This result is not an approximation. There are already nu- 
merous statements in the literature that the actual gravitational 
field may in some cases be approximated by a uniform gravita- 
tional field without introducing large errors. Such conclusions 
might be based, for example, on series expansion procedures, and 
indicate that the perturbation problem is not singular. In con- 
trast to these statements the conclusions given in this report are 
intended to be exact. 

The results are not restricted to planar problems, and ex- 
tensions to non-central and time-dependent fields can be made 
without difficulty. 

A material increase in the information obtained would result 
from the use of better bounds on the useful velocity perturbations 
which can be produced by the force field. Further investigation 
of these bounds might be profitable. 


(2) Introduction 


ge involving the transfer of a rocket vehicle 
from one point to another in a gravitational field 
with minimum fuel expenditure have been studied by 
Lawden,! Fried,? and many other authors. Generally, 
the methods of variational calculus have been employed, 
and many valuable results have been obtained. How- 
ever, it seems that the subject is sufficiently complex 
* to warrant attempting other procedures which may also 
: contribute to the understanding of these problems. 

2 In the present report a known condition for optimum 
transfer in uniform gravitational fields (without at- 
mospheric resistance), that the impulses be applied 
only at the terminals,* is extended to a class of problems 
in centrally directed fields. 


Received March 13, 1959. Revised and received July 28, 1959. 
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This extension is not an approximation, but is in- 
tended to be exact. There are already numerous state- 
ments in the literature (see references 4 and 5, for ex- 
ample) regarding the approximation of nonuniform 
gravitational fields by uniform fields. Such results 
might arise in series expansion methods, for example, 
and indicate that the perturbation problem is not a 
singular one. The present developments are exact 
even though a linearized force field is employed, be- 
cause the optimum trajectory is compared only with 
infinitesimally different paths. This will be discussed 
in more detail later. 

The known results for force-free and uniform gravi- 
tational fields are first redeveloped by using an acceler- 
ating but nonrotating coordinate system whose origin 
follows a coasting trajectory in the given gravitational 
field. The same method is then applied in centrally 
directed gravitational fields. 

For clarity of presentation most of the development 
is given for “‘two-dimensional”’ cases, in which the coast- 
ing trajectory and terminal velocity vectors are all 
contained in one plane. However, the conclusions 
are readily extended to “‘three-dimensional” problems. 


(3) The Force-Free Space 
(a) Statement of the Problem 


Consider the following problem: given a particle 
traveling in a two-dimensional force-free space. At 
time ¢ = 0 the particle is at point Py’ (see Fig. 1) 


y’ 


Fic. 1. A typical problem in the fixed coordinate system. 
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Fic. 2. A typical problem in the moving coordinate system. 


traveling with vector velocity Vo’. At time ¢ = 7 the 
particle must be at point Pr’ traveling with vector 
velocity Vy’. What is the procedure for accomplishing 
this with minimum impulse? 

It can be shown that the minimum impulse is at- 
tained when the sum of the magnitudes of the velocity 
changes is minimized. (This is true even when the mass 
is decreasing due to fuel expended. See, for example, 
reference 1.) We then wish to minimize AV where 


AV = |AV,| + |AV.| + +... (1) 


Each AV is assumed to be impulsively generated at a 
particular time (no two occurring at the same time). 

It is useful to reformulate this problem in a non- 
rotating but uniformly translating coordinate system. 
Let the moving coordinate system have velocity * = 
d'/T relative to the fixed coordinate system and move 
parallel to the straight line connecting Po’ and Pp’ 
and in the direction Py)’ ~ Pr’. In the moving co- 
ordinate system (see Fig. 2) the initial and final points 
then coincide. The terminal velocity conditions have 
been altered by the velocity of the moving system, but 
a particle still encounters only such forces as are arti- 
ficially applied (the space is still force-free). 

The problem may now be stated as follows: at time 
t = 0a particle is at the origin traveling with vector 
velocity Vo. At time ¢ = 7 it must again be at the 
origin but traveling with vector velocity Vr. How is 
this to be accomplished with the minimum AJ’? 


(b) Single Intermediate Impulses 


We now consider the possibility of applying velocity 
changes AV;, AV; at times ¢ = 0, ¢ = 7, and t = 
t;, respectively, where ¢; is an arbitrary intermediate 
time. If u denotes the x component of velocity and v 
the y component of velocity in the moving coordinate 
system, then a typical velocity diagram in the u, v or 
hodograph plane will be as shown in Fig. 3, at the left. 
The corresponding path of the particle in the x, y plane 
is shown at the right. 

It is important to note that the velocity vector AV, 
must always pass through the origin in the u, v plane. 


The particle leaves the origin of the x, y plane in a di- 
rection shown by the vector sum of Vy and AVo, and 
travels in a straight line with constant speed. The 
statement of the problem requires that the particle 
must return to the origin of the x, y plane by time ¢ = 
7. Thus, AV; must have a direction opposite to that 
of Vo + AV» and an equal or greater magnitude in order 
to cancel the original outward velocity and create an 
inward motion. 

The total AV required for any chosen path from A 
to B in the u, v plane is equal to the length of that path 
(a scalar quantity). The minimum path length from 
A to B passing through the origin is |Vo! + |Vr|. In 
general this can only be attained by setting AV) = 
—Vo, AV; = 0, and AV; = Vy. In the moving coordi- 
nate system the particle then stays at the origin from 
t = Ountilt = 7. This means that the optimum pro- 
cedure is to use impulses at the terminal points with a 
coasting trajectory between terminals. 

However, when V) and Vy are directly opposed 
the minimum value of AV can also be attained in 
many other ways with AV, ¥ 0. 

So far we have considered only a “two-dimensional” 
problem. In a “three-dimensional’’ case we will let 
the x, y plane of the moving coordinate system contain 
the coasting trajectory between the terminal points. 
The terminal velocity vectors may then have com- 
ponents normal to the x, y plane. However, the short- 
est path between two points in the hodograph space 
(u, v, w) by way of the origin is still composed of two 
straight lines. 

The following statement has then been proved: 
given a particle traveling in a force-free three-dimen- 
sional space. At time ¢ = 0 the particle is at point 
P,' traveling with vector velocity Vo’. At time ¢ = T 
the particle must be at point Py’ traveling with vector 
velocity Vr’. Consider the possibility of applying three 
impulsive thrusts, one at ¢ = 0, one at ¢ = 7, and one 
at ¢ = ¢, an arbitrary intermediate time. Then 
eliminating the intermediate impulse never increases 
the minimum total impulse required. 
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Fic. 3. A typical excursion from the origin in the moving co- 
ordinate system. 
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Fic, 4. A path in the fixed coordinate system with multiple 
intermediate impulses. 


(c) Multiple Intermediate Impulses 


So far we have considered only the possibility of 
impulses at = 0,¢ = 7, and?t = t;. Suppose we now 
permit impulses at any finite number of intermediate 
times, so that in the original fixed coordinates the path 
between P)’ and Pr’ might consist of many straight 
line elements as shown in Fig. 4. 

We take any three consecutive points, k, k + 1, and 
k + 2, and assume the terminal vector velocities speci- 
fied at k and k + 2 and the time of transit fixed. Our 
previous conclusions indicate that the same or less 
impulse is required in going directly from k to k + 2 
without the intermediate point k + 1. Repeated 
application of this method allows us to eliminate all 
intermediate points between Po’ and Pr’ and conclude 
that the optimum path is a straight line with impulsive 
thrusts applied only at Po’ and Pr’. The illustration, 
Fig. 4, shows a two-dimensional path, but the principle 
applies equally well in three dimensions. 

A similar result can be obtained by the use of super- 
position methods, which are discussed later in connec- 
tion with the centrally directed gravitational field. 


(4) The Uniform Gravitational Field 


If we have a uniform gravitational field instead of a 
force-free field, a translating (nonrotating) coordinate 
system with constant acceleration can be used. First, 
in the fixed coordinate system we have points P)’ 
and P,’ as shown in Fig. 5. 

The uniform gravitational field is such as to produce 
an acceleration ‘‘g’’ in the negative y’ direction. 
Again we prescribe a vector velocity Vo’ at Po’ at t = 0 
and a vector velocity Vr’ at Pr’ att = 7. Let the 
moving coordinate system have a constant horizontal 
velocity V,*. 


V,* = (xp’ — x0’)/T (2) 


and a variable vertical velocity V,* 
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V,* = [(yr’ — yo’ + gT?/2)/T] — gt (3) 


Then in the moving system the initial and final points 
will coincide. Also the moving coordinate system is 
force-free since the gravitational force is balanced by 
the effective force created by the constant acceleration. 
The problem in this moving system is then exactly 
the same as that studied previously, and the same con- 
clusions follow. 


(5) The Centrally Directed Gravitational Field 


(a) The Measure of Nonuniformity 


Let uw define the strength of a central force field 
through the relation 


central force per unit mass = 
(radius from center)? 


Let 7’ be the radius from the center of force to the 
origin of the translating coordinate system, and let 
T be the time specified for the transfer. Then the 
dimensionless parameter P, defined by 


Py = (5) 


measures the importance of the local nonuniformity of 
the central force field. This parameter will appear 
naturally in later developments, where its most critical 
value, P, will be of interest. It measures the velocity 
increments which can be produced by the perturbation 
force field. (This force field, which is effective in the 
moving coordinate system, is closely related to the force 
field which produces tides on the earth.) 

If uw = Qorr’ = ~ or J = 0, then we have P,; = 0. 
The gravitational field is then uniform, and the per- 
turbation force field, being identically zero, produces 
no velocity increments. 

In general, since the origin of our moving (but non- 
rotating) coordinate system follows a coasting tra- 
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Fic. 5. A problem in the uniform gravitational field. 
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Fic. 6. Coordinate systems used in expressing the perturbation 
force field. 


jectory with angular velocity w’ about the center of 
force, 


In the special case when this trajectory is an are of a 
circular orbit, 7)’ = O and 


= w"? (7) 
then in this special case 
P, = oT? = @” (8) 


where 0’ is the angle corresponding to the arc traversed. 

We will be concerned primarily with cases in which 
P, < 1.0. Note that this does not necessarily restrict 
the results to small angular displacements about the 
center of force. If the specified time for the transfer 
T is small, then the results may apply even when large 
angular displacements (of order 180°) about the center 
of force occur. It can be seen that when 7° is very 
small the impulses contributed by gravitational forces 
during the transfer (being proportional to 7) may be 
negligible compared to the other impulses involved. 
This can be true even when the angular displacement is 
large. In such a case the coasting trajectory would not 
be a circular are about the center of force and the rela- 
tion P, = 6” would not apply. 


(b) Single Intermediate Impulses 


We will use a nonrotating coordinate system whose 
origin follows a coasting trajectory in the given gravi- 
tational field, as shown in Fig. 6. 

We wish to find the maximum magnitude of force 
which can exist at radius 7 from the origin of the moving 
system. At the origin the force must, of course, be 
zero since the origin follows a coasting trajectory. 

First we determine the x and y components of the 
gravitational force field. 


(9) 


F, = —(mp/r’?) cos 
F, —(mp/r’?) sin 


At the origin of the moving coordinate system these 
forces are exactly balanced by the accelerations ex- 
perienced by the coordinate system. Since the coordi- 
nate system has only translational motions these 
accelerations are independent of position in the moving 
coordinate system. However, the gravitational force 
field varies with position, and this creates the perturba- 
tion force field. We plan to point out the situations 
in which only terminal impulses are desirable. As 
in the calculus of variations we need compare only 
paths lying in the immediate neighborhood of the opti- 
mum, which is the coasting trajectory. Thus, all ex- 
cursions will remain in the neighborhood of the origin 
of the moving coordinate system, and we can use the 
linearized perturbation force field. For this we need 
the partial derivatives of the force components with 
respect to displacement evaluated atx = y= 0. These 
are as follows: 


OF,/dx = —my}(1 — 3 cos? 

OF,/Oy = —my}(—3 sin cos A") /ro’*} 

OF,/Ox = —my}(—3 sin cos 6’) ro'3} (10) 
oF, = —my} (1 — 3 sin? &’) /ro’*t 


Here 7’ and @’ are functions of time determined by 
the coasting trajectory which is followed by the origin 
of the coordinate system. 

The motion of a particle in the linearized force field 
is then defined by 


(1 — 3 cos? 0’) — 
uy| (—3 sin cos 6’) 
—px}(—3 sin cos — 
uy} (1 — 3 sin? 


(11) 


Xe: 
Il 


Let x = r cos 6, y = r sin 0, then 


(ur/ro’*)} [3 cos? 6’ — 1] cos + 


[3 sin 6)’ cos | sin 6} (12) 
= (ur/ro’)} [3 sin? — 1] sin 


[3 sin @)’ cos cos 6} 


We wish to find the maximum possible acceleration 
which can be obtained at a given time. To facilitate 
this we choose our system of axes so that #’ = 0 at that 
time. Then we get 


= (ur/ro’*)2 cos 6} | 
¥ = (ur/ro’’)} —sin 6} 


(13) 


The maximum resultant acceleration is then 


Amar = Qyr/19"* (14) 


This was derived in two dimensions, but applies 
equally well in three. The radius r is then the radius 
in x, y, Z space, not necessarily lying in the x, y plane. 

It is now necessary to find the magnitude of the per- 
turbation velocities which can be produced in any ex- 
cursion from the origin of the moving coordinate system. 

For the outbound trip we define a mean vector ve- 
locity V;, and similarly for the inbound trip V2,,, such 
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that the following relation applies: 
Vim 6h, = R = — Vo», bto (15) 


where R is the position at the time when the intermedi- 
ate impulse is applied. Also, 


8t, + bf = T (16) 


If the radius R at the time when the intermediate 
impulse is applied is also the maximum radius attained 
(Rar) then an upper bound on the magnitude of the 
acceleration is 


a= (Qu ate (17) 


We will assume temporarily that R is Ry», then see 
under what conditions this is true. It is desirable to 
find an upper bound on the magnitude of the deviation 
of any instantaneous velocity from the mean velocity. 
This should occur when the incremental velocity vectors 
all have the same direction, and the upper bound on 
the acceleration is applied at all times. For this linear 
variation of velocity components with time the mean 
velocity is that occurring at the half-way point in time. 
The greatest departure from the mean velocity on the 
outbound trip is then not greater than 


6V, = (1/2)a bt, (18) 
6Vi/Vim = (u/ro min’®)6t,? (19) 

and on the inbound trip 
5V2/ Vom = min’*) (20) 


We introduce here and call it P. It can 
be shown that R is Kya: if P < 1.0 by studying a suc- 
cession of points along a hypothetical path in the x, y 
plane. We start at a point on the path where the local 
radius 7, reached at time ¢, is the maximum so far 
attained in moving out from the origin. Then the 
6V accumulated at 7 (measured from the mean ve- 
locity r/t) is bounded by 


8V < (1/2)at (21) 
Since @ = 2ur/romin’® we get 


< urt/ro (22) 
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Fic. 7. Path in the physical plane and restriction on velocities 
for F<. 1.0. 
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Fic. 8. Comparison of path lengths in the hodograph plane. 


P = pT?/ro min’? and P will be assumed < 1.0. Then 
eliminating min’? we have 


6V <rt/T? or 6V/(r/t) < /T? (23) 
and éV/(r/t) < 1.0 (24) 


This means that the instantaneous velocity at time 
t always has a component directed radially outwards, 
as can be seen in Fig. 7. A similar analysis can be 
made for the inbound portion of the trip. It is then 
clear, through application of this argument to all 
points on the path, that RX must be R,,,, under the re- 
striction that P < 1.0. It follows that the expressions 
derived for 6V1/Vim and 6V2/V2m are correct for P < 
1.0. 

The bounds just obtained on 6V; and 6V2 are not the 
best possible, and it may be desirable to improve them 
so that more information can be obtained from the 
final equations. For example, in the present develop- 


ment it has been assumed that the maximum accelera- | 
tion which can be developed at the extreme radial dis- | 


tance & acts at all times, although the average accelera- 
tion must be appreciably less than this. Also, it has 


been assumed that the maximum acceleration is avail- | 


able in any desired direction. Actually, the accelera- 
tions available have specific directions ranging from 
inward to outward, and the latter have a maximum 
twice as great as the former. At the expense of addi- 


tional complexity in the analysis, corresponding to | 


more accurate evaluation of these factors, the bounds 


on 6V; and 6V2 might be improved. However, this © 


will not be attempted at present. 

In order to determine when an intermediate impulse 
cannot be used to decrease the total impulse, we com- 
pare the paths indicated by the solid line and the 
dashed line in the hodograph plane, Fig. 8. The dashed 
line represents the use of terminal impulses only. The 
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solid line corresponds to the introduction of an inter- 
mediate impulse, and is chosen as the shortest line 
which connects A with B by way of the two circles. 
The vector AV, then passes through the origin, since 
Vim and V2, are opposite in direction. (In a specific 
example, using the equations of motion instead of 
bounds on the velocity increments, AV; would not 
generally pass through the origin.) 

The upper figure is essentially an enlargement of the 
lower one. The perturbation from a single coasting 
trajectory is assumed small so that the outer portions 
of the solid and dashed lines become parallel. This is 
consistent with the fact that a linearized perturbation 
force field is used. Thus, as in the calculus of vari- 
ations, we are comparing the optimum path with neigh- 
boring paths, not widely different ones. 

The circles indicate upper bounds to the velocity 
contributions of the perturbation force field. The 


path length for the solid line cannot be shortened more 
than 2(6V, + 6V2) from the length of the three con- 
secutive straight line segments on which the solid line 
is superimposed. It is obvious that if 6V, and 6), 
shrink to zero the solid line can never be shorter than 
the dashed line. If 6V, and 62 are small, but not 
zero, there will be a large range of angles y for which 
the solid line cannot be shorter. 

The parameters 8, V;,,, and V2,, determine the direc- 
tion and time of application of the intermediate impulse. 
We wish to use the solid line (or intermediate impulse) 
in the best possible fashion, so must adjust @ and 
V,, and V2,, to their optimum values. We then find 
the conditions under which the two path lengths are 
exactly equal. This should define some sort of bound- 
ary, on one side of which the use of terminal impulses 
alone will always be optimum. 

Neglecting the parallel portions of the paths, which 
are of equal length, we equate the path lengths and get 


Vim — 6Vi + Vom — V2 = Vim sin [y/2 — B] + EV; + Vom sin [y/2 + B] + 6V2 (25) 


Replace 6V, by its upper bound, min’®) and 6V2 by its upper bound, in’), 


Vim {1 — ‘ro — sin (y/2 — B)] + Vom{1 — (2u/ro min’®)dte2 — sin (y/2 + B)| = 0 (26) 

We also have the relations 
Vom = Vim 6t,/dte (27) 
and 64, + bt2 = T (28) 


Eliminating V2,, and 6t, by means of these relations gives 


(T — (1 — (24/70 — sin (y/2 — B)] + — (2u/r0 mn’®)(T — — sin (y/2 + B)] = 0 (29) 


Let 6t,/T = (1 + £)/2 (30) 
Then T — & = T(1 — £)/2 (31) 

(1 — &)[1 — (2¢7?/ro min’®)(1 + &)?/4 — sin — + 
(1 + &)(1 — (2u7?/ro min’®)(1 — &)?/4 — sin (y/2 + B)] = 0 (32 
ul” = 2 — (1 — &) sin (y/2 — B) — (1 + &) sin (7/2 + 8) (33) 

Expanding sine terms and replacing u7™*/ro min’* by P gives 

p=2 — sin (y/2) A cos /2) (34) 


This shows the origin of the parameter P which meas- 
ures the nonuniformity of the gravitational field. 

At this point we have a relation between P and y 
which depends on & and 8. When this relation is 
satisfied the use of terminal impulses alone and the use 
of terminal impulses plus an intermediate impulse re- 
quire equal total impulses under our assumptions. 
We wish to be able to say that for a given value of y if 
P lies below a certain value then terminal impulses 
alone are always better. This means that for a given 
value of y we must find the minimum value of P which 
is attainable by varying & and 8. 

If we differentiate P with respect to 8 and set the 
derivative equal to zero we get 


£ = tan B tan (7/2) (35) 


while differentiating with respect to £ yields 


(1 + &) sin B cos (7/2) = 
2¢[1 — cos B sin (y/2)] (26) 


Eliminating either or 8 and solving, we find that 
the stationary values of P occur for 


B=0, €=0 
B= (nr —y)/2, €=1.0 (37) 
B= — y)/2, —1.0) 


The minimum value of P occurs when 8 = +(a — y)/2 
and = +1.0 and is given by 


P = cos? (7/2) (38) 
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Fic. 9. The P versus y diagram. (vy is the lesser angle between 
the two vectors measured in their common plane.) 


(Since 6t, lies between 0 and 7, — lies between —1 and 
+1. If stationary values of P had not appeared at 
& = +1 these end points would have required separate 
investigation to see if the minimum value of P occurred 
there.) 

We now plot P versus y in Fig. 9, and indicate, by 
shading, the area which yields useful information. 
The boundary corresponds to 


Prin = cos? (y/2) 


At the moment we have shown that terminal impulses 
are optimum in the shaded area for two-dimensional 
problems in which only single intermediate impulses 
are considered. Extensions to three-dimensional prob- 
lems and multiple intermediate impulses are discussed 
below. 


(c) Extension to Three Dimensions 


So far we have only compared path lengths in the u, 
v plane. We can also consider cases in which the 
vectors V;,, and V2, have components normal to the 
u,v plane. In Fig. 8 it can be seen that the length of 
the solid line must increase in such cases. Thus, we 
have already made the best possible use of the solid 
line (or the intermediate impulse). 

The general extension to three-dimensional problems 
follows from this. In such problems the vector veloci- 
ties prescribed at the terminal points in space do not 
necessarily lie in the plane of the coasting trajectory 
connecting the points. We may then have Vp and 
Vr not contained in the u, v plane. However Vp and 
Vr always define some plane in the hodograph (uw, v, w) 
space, and in that plane the previous developments 
apply exactly as before. 


(d) Multiple Intermediate Impulses 


The equations of motion governing an excursion 
in the x, y plane (the moving coordinate system) have 
been linearized. This restricts our study to excursions 
which remain in the neighborhood of the origin, as 
previously mentioned. However, the linearity makes 
superposition possible. The form of the equations 
[see Eqs. (11)] is 
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= falt)x + + (39) 
= + + Ba(t,)! 


where x and y are functions of time ¢, # denotes d?x/dt?, 
etc., and 4(t;) is the Dirac delta function corresponding 
to the intermediate impulse. 

Let one excursion in the x, y plane be given by x,)(t), 
y¥ay(t) and a second excursion by Y(t). Then 
Xay(t) + x(t), + satisfies the equations 
of motion. If the intermediate impulse was applied 
at time in Xq)(4), Ya) (t) and at time in (0) 
then the combined excursion will contain two inter- 
mediate impulses. In the same fashion multiple inter- 
mediate impulses can be obtained, and all multiple 
intermediate impulse cases can be decomposed into 
single impulse cases. 

It follows also from linearity that path lengths in the 
hodograph plane are additive for two or more separate 
excursions, and that the reduction in path length which 
might be obtained with multiple intermediate impulses 
is the sum of the reductions in path length obtainable 
with the component single intermediate impulse ex- 
cursions. 

From this we conclude that if there is no way to 
apply any single intermediate impulse and reduce the 
total impulse required, then multiple impulses are 
also ineffective. This generalizes our results. (The 
same approach can, of course, be applied in the force- 
free field and uniform gravitational field.) 


(6) Results and Conclusions 


(a) The diagram in Fig. 9 may now be used as 
follows: 

Given two points Po’ and P,’ in a central force field, 
specified velocity vectors Vo’ and Vr’ (any direction) 
at the points, and a specified time, 7, for the transfer 
between points. 

First, find a coasting trajectory between the two 
points which requires time, 7, for transfer. (Generally 
there are two such trajectories, which pass on opposite 
sides of the center of force. Frequently only one of 
these will be of interest.) 

Second, find the minimum radius, 7 min’, for this 
trajectory. 

Third, compute P, where P = min’. 

Fourth, find the terminal velocity vectors Vp and 
Vr in the moving coordinate system (a nonrotating system 
whose origin follows the coasting trajectory). Vp is 
then the vector difference Vo’ — Vj,aj, at Po’ and Vr is 
Ve’ — Viraj. at Pr’, where V;,;, is the vector velocity of 
a particle following the coasting trajectory. 

Fifth, find the (smaller) angle, y, between Vp and Vr 
measured in their common plane. 

If the point defined by P, vy falls inside the shaded 
area of Fig. 9, then the minimum possible total impulse 
for transfer (subject to the terminal velocity condi- 
tions) is obtained by using impulses only at the termi- 
nals. The use of intermediate impulses may not in- 
crease the total required impulse (in special cases) but 
can never reduce the total required impulse. 
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If the point, P, y falls outside the shaded area then no 
positive conclusion can be drawn. A more detailed in- 
vestigation must be made to determine the optimum 
procedure. 

(b) Fig. 9 applies only to a central force field. _How- 
ever, similar methods can be developed to cover more 
complicated fields, including those which are time 
dependent. It is then necessary to find an upper bound 
on the acceleration which can be encountered in the 
moving coordinate system for such cases. 

(c) It seems probable that the shaded area in the 
plot of P versus y (Fig. 9) can be enlarged by obtaining 
better upper bounds to the useful velocity increments 
which can be contributed by the perturbation force 
field. This would increase the information obtainable, 
and is worth investigating. 

(d) Finally, we emphasize that the results obtained 
are intended to be precise. The analysis employs upper 
bounds on the effects of the perturbation force field, 
not merely approximations. This limits the amount 


of information obtained, but those results which do 
appear should be exact. The linearization of the 
perturbation force field is justified by the fact that 
we study only infinitesimal excursions from the origin. 
The optimum trajectory (in this case the coasting 
trajectory) is then compared only with neighboring 
paths, not widely different ones. 
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powerful analytic tools of hydrodynamics can be made 
directly applicable to the shear lag problem. The 
analogy also suggests the use of graphical network 
solutions when highly unusual boundary configurations 
must be dealt with. It is also conceivable that at cer- 
tain times it might be advantageous to solve shear lag 
problems by direct experimental observation of the 
corresponding plane problem of fluid flow. 

Finally, in addition to its value for numerical prob- 
lem solving, it is believed that the analogy affords the 


designer a simple physical picture of the stress flow 
which can be used to advantage for intuitive appraisals 
of structural integrity. 
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Natural Frequency of Beams With Variable 
Moments of Inertia 


Ming L. Pei 

Associate Professor, Civil Engineering Department, 
City College of New York, N.Y. 

November 20, 1959 


= DIFFERENTIAL equation governing the free flexural vibra- 
tion of an elastic beam with variable flexural rigidities E/(x) 


and nonuniform mass per unit length m(x) 


(Ely")” = w*my 


(1) 


may be expanded into finite-difference equations and expressed 


in matrix notations as! 


[B] [EI] |B) = 


(2) 


where 

[Y] = column matrix of deflections 

[EI] = diagonal matrix of flexural rigidities of the beam at 
n — 1 equally spaced points 

[M] = diagonal matrix of mass per unit length at n — 1 
equally spaced points 

h = L/n, the distance between two adjacent points 

w = circular frequency 

{B] = finite-difference operator for d?/dx? 


In the case of simply supported beam, yo = 
isan (m — 1) X (m -- 1) square matrix as follows: 


1 1 
1 -2 1 
= 1/h 1 
1 
1 


yn = 0; and [B] 
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[B] 


(3) 


Though Eq. (2) can be solved by iteration process, unfortunately, 


the highest, rather than the fundamental, frequency is ob- 
tained.?, In order to obtain the fundamental frequency as the 
first result of the computation, it is necessary to modify Eq. (2) 
to the following: 


[Y¥] = w?[B] ~'[h/EI) [M][Y] (4) 
or, simply, 
[Y] = w?[C] [Y] (5) 
where 
= [h/EI)|B)— (6) 
[h/EI) = = 


a diagonal matrix with h/EJ; as diagonal elements 


The fundamental frequency is obtained when Eq. (5) is iterated. 
However, the scheme is little used because of the computational 
efforts required for the inversion of [B], particularly when n is 


large. It is noted that |C] is the flexibility matrix of the beam 
and is sometimes obtained directly by a static analysis. This 
requires a similarly lengthy calculation when n is large. It has 


been pointed out! that, surprisingly, the inverse of [B] can be 
expressed by a formula of closed form. For simply supported 
beam, 
[A] = [B]™ (7) 
(i(n — j)h/n i<j 


lin (8) 


ay; = 
For any given distribution of mass and EJ, we can compute 
|C] with the aid of Eq. (8). Knowing [C], iteration of Eq. (5) 
will yield the desired fundamental frequency and mode of the 
beam. The method does not require matrix inversion and is 
very suitable for digital computers, as the programing contains 
nothing else besides a few basic matrix operations. 
The method can be used for cantilever beam also. With the 
beam fixed at x = O and free at x = nh = L, Eq. (4) becomes, 


[Y] = w?[D]’[h/EI) (9) 
where 


[D]’ = transpose of |D] 
[D] = (n+ 1) X (n + 1) square matrix 


= (G- ah i>j; j<n (10) 
— 1)h/2 i>j; j=n 
0 i<j 


It should be noted that the inversion of difference operator [B] 
is equivalent to a double integration process. The fact that the 
resulting matrices [A] and |D] are different is due to different 
boundary conditions which have been taken into account already. 
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Slip Coefficient for a Tire 


John E. Stevens 

Chief of Space Mechanics, Astronautics Division, 
Chance Vought Aircraft, Inc., Dallas, Texas 
November 6, 1959 


(1) INTRODUCTION 
T IS WELL KNOWN that when a brake torque is applied to a 
wheel the rotational velocity of the wheel slows down to a 
value less than would be expected from the forward speed of the 
wheel and the rolling radius of the wheel. This phenomenon is 
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explained by the strains set up in the tire wall due to the drag 
force. The strain is proportional to (Aw) R/V, where R is the 
rolling radius, V is the forward velocity, and (Aw) is the change 
in rotational velocity. If the friction coefficient is low, it is 
assumed here that the slip situation may be explained by assuming 
that the part of tire in contact with the ground is skidding at all 
times and the tire walls are effectively rigid. This development 
makes use of the fact that the elements of the tire which come 
in contact with the ground each has a different rolling radius. 
The elements of the tire along the centerline of the tire are going 
faster than those nearest to the edge of the ground contact area. 
Since the net drag force without braking is nearly zero, this 
means that the elements are going the fastest slip aft giving rise 
to a forward drag force and those going the slowest slip forward 
giving rise to an aft drag force. When a drag force is applied 
to the wheel, the number of elements slipping forward increases. 
The number of elements slipping forward compared with the 
number slipping aft determines the average drag force and the 
average slip velocity. 


(2) THEORETICAL DEVELOPMENT 
It will be assumed that the average friction force on each 
differential area of tire in contact with the ground can be written 
as 
(dF/dA)p = wP,/rOprho 
where (see Fig. 1) 


the friction coefficient 


= 
P, = the vertical load on the tire 
Orr = half the footprint width 
hy = half the footprint length 


The net drag force on the tire is, then, 
1.0 


V1 — (6/0)? — 


| V1 —(0/6R)? aon) | 
0 


Fp = 


and 
dFp/d0 = — (6/0n)? 


where @ is the angle where the slip force changes sign. 
velocity change with @ can be written as 


dv/d@ = —w(d/dé)[R — r(1 — cos 6)] 


The slip 


= wr sin 6 


~ 
| 
(-} 
Ta- | 
(x) | 
(1) | 
sed 
(2) 
at 
1 | 
j 
& 


dv 


_ _ V1 — (6/0r)? 
do dv rVrOp sin 6 


where V is the forward speed. 


6 can be determined by setting Fp = 0 and solving for @. 
= sin~! (6/0r) + 1 — (0/Or)? 


or 
= 
From reference 1, 
~ 2 »/(8/w) — (8/w)? 
where W is the width of the tire and 6 is the vertical deflection of 
the tire. 

If we define K, by 

Fp = K.AV/V = K.(dwR/V) 
we can then write 

K, ~ 16DpPv/wr?[6/w — (6/w)?] ~ 1.6DyuP./5 = 1.6DuK, 
where K, is the vertical spring rate of the tire and D is the di- 
ameter of the tire. 

If the test data on the slip behavior of tires on dry concrete 
from reference 1 are compared with the value of the slip coeffi- 
cient developed here, it will be seen that, if dry concrete values 
for the sliding friction coefficient are used, this development pre- 
dicts slip coefficients significantly greater than those measured. 
However, if the surface condition is such that the friction co- 
efficient is low enough, the mechanism described here may offer 
a better description of the physical situation. 
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Experimental Study of an Ablating Sphere 
With Hydromagnetic Effect Included} 


John H. Boynton* 
Research Engineer, Convair-Astronautics, San Diego, Calif. 
November 12, 1959 


SYMBOLS 
B = induced magnetic field, scalar 
B = potential-velocity constant 
Cw = specific heat of water 
im = thickness of the liquid-metal layer 
H = applied magnetic field, scalar 
Lan = heat of fusion of the alloy 
pw/pm = density ratio, water to metal 
Ss = conductivity of the alloy, electrical 
ow = Prandtl number of water at the wall 
To, Tw = temperature, free-stream and at the wall 
t = time 
a = average velocity in the liquid-metal layer 
um = absolute viscosity of the liquid metal 
vw, ¥m = kinematic viscosity, water and liquid metal 
Wm = melt rate of the alloy 
Zs = stagnation-point position 


NTEREST in solving the aerodynamic heating problem associated 
with re-entry into the earth’s atmosphere has led scientists to 
develop the mechanism of material ablation and, more recently, 
to study theoretically magnetohydrodynamic interaction with the 
ionized shock layer. The present study! was conducted in an 
effort to verify theories explaining each of these effects, and the 


+ This research has been supported in part by the OSR of the USAF under 
Contract No. 18(600)-1523. 
* Formerly, Research Assistant, Cornell University. 
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possibility of a re-entry vehicle employing both of these phe- 
nomena is offered. Ablation as an energy-dissipation device has 
been operationally proved, and theoreticians?~* have shown that 
heat-transfer reductions of 30 per cent may be realized under 
normal re-entry conditions for the case of an applied magnetic 
field. 

The hemispherical models were cast from an eutectic alloy of 
bismuth, and the low melting temperature of 117°F. permitted 
the use of hot water as the flow medium. The apparatus con- 
sisted essentially of a ‘“‘blow-down’’ water tunnel with a 5-in. 
vertical test section. Data were recorded on 16-mm. film, as 
depicted in Fig. 1, and analyzed through the use of an optical 
comparator. Instrumentation included a stopwatch, free-stream 
thermocouple, and a water-level gauge. For the magnetic portion 
of the experiment, a coil around the model support in conjunction 
with a six-volt storage battery yielded a field estimated at 5,400 
Gauss. The relatively high electrical conductivity of the alloy 
(4 per cent of Cu) permitted a hydromagnetic interaction be- 
tween the applied field and the liquid-metal layer. 

Assuming potential flow outside the boundary layer, Sibulkin® 
has derived an expression for the heat transfer at the stagnation 
point of a sphere, with the analysis based on a Nusselt number 
relation obtained by numerical integration. Balancing the heat 
input to the body with that used to produce the melted layer, a 
nondimensional melt rate’ may be formulated as a function of 
free-stream temperature: 


The experimental ablation results for various free-stream temper- 
atures are plotted, together with Eq. (1), in Fig. 2. 

The magnetic circuit was constructed such that the only air- 
space was between the tip of the model support and a ring around 
the test section, thereby offering the most efficient geometry anda 
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Fic. 1. Sequence taken from 16-mm. film of ablation experi- 
ment. Time increases from top left to bottom right as indicated 
by stop watch. Note that model retains smooth, approximately 
hemispherical shape during ablation but diminishes in size. A 
drop of molten metal can be seen leaving the model skirt in the 


fourth frame of the sequence (middle right). 
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Fic. 2 (top). Melt rate vs. temperature difference. 
Fic. 3 (bottom). Effect of magnetic field on ablation. 


field that was nearly radial at the body surface. An equivalent 
rectangular airgap was used to estimate the field strength at the 
model, and, since the pertinent magnetic interaction parameter, 
H5m/umt, was greater than unity, a discernible magnetic effect 
could be predicted. This interaction essentially increased the 
effective viscosity of the molten layer, according to the relation 


Vett. = + (Siim )( B75 (2) 


which retarded the flow and thickened the liquid-metal layer. 
This appears to be a decrease in the melt rate in Fig. 3, but, since 
the heat-transfer process depends almost wholly on conduction, 
it is actually a thickening of the molten layer. The quantity 
plotted against time in Fig. 3 is the stagnation-point position 
from a fixed reference point. 

The melt rate without the magnetic field is given, within ex- 
perimental error, by Sibulkin’s incompressible heat-transfer 
analysis. A qualitative effect is observed when the magnetic 
field is applied; the liquid layer appears to thicken. This result 
indicates the possibility of efficiently reducing heat-transfer rates 
to a blunt body by increasing the amount of material which 
vaporizes in re-entry. The results also indicate the need for 
further quantitative experiments with possible comparison with 
appropriate theory. 
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The Transfer Function of a Spring-Tab 
Control 


Malcolm J. Abzug 
Douglas Aircraft Company, Inc., El Segundo, Calif. 
November 6, 1959 


SYMBOLS 
Ch = H/Stq 
Chg = OCh/ 05 
= 0Cp/0(d5/dt) 
Chg, = OCh/ O5¢ 


(Chg) yao = (0H/06, tab locked, at zero airspeed) /Stg 
(Chg, peo = (0H/0é:, control surface locked, at zero airspeed) /Siq 


Ci = I/Stq 

H = control-surface hinge moment 

K = gain 

= dynamic pressure 

Ss = control-surface area 

s = Laplace-transform variable 

t = control-surface chord 

6 = control-surface angle 

bt = tab angle 

5t5 = 065;/05, input locked 

= control surface locked 

t = damping ratio 

o = input (auto-pilot servo or pilot’s cockpit control) angular 
deflection 

@n = undamped natural frequency 

(wn)a = undamped natural frequency of airplane’s short-period mode 


of motion 


A POINTED OUT by Passera and Nason! and by White,? the 
performance of closed-loop flight-control systems is often 
critically dependent upon the dynamics of the servos used for 
control-surface actuation. The dynamics of two hydraulic 
servomotors typical of control-surface servos have been discussed 
by Gold, Otto, and Ransom.’ It is the purpose of this note to 
examine the dynamics of a well-known type of aerodynamic 
servo, the spring-tab control, from the standpoint of closed-loop 
flight control. A recent paper by Adamson and Lyons‘ has 
treated the dynamics of the spring-tab control as they affect 
manual control in several special maneuvers. 

The control under consideration is the typical geared or un- 
geared spring-tab control described by Phillips? and Dunn.* 
The input to the spring-tab control is obtained from either an 
auto-pilot servo or from the pilot’s cockpit controls. The spring- 
tab and control-surface combination forms one of an array of 
dynamic elements, a block in the closed-loop block diagram. An 
equation of motion for the control surface about its hinge line 
and a kinematical equation describe the system: 


[Cig + + [Cng, + (Crs, + 
Cys*6 = 0 (1) 
br = + (2) 


The overall spring-tab control transfer function is evidently the 
output control-surface deflection 6 divided by the input deflection 
oa. Solving Eqs. (1) and (2) simultaneously: 
6 1 
= K (3) 
o (s/w,)? + 2fs/w, +1 
where 


ig [Cy + + 


= 
| 
peri- 
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he 
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TABLE 1 
‘ 
Transfer Function Properties of Typical Spring-Tab Controls 
FLIGHT PHASE LAG 
1GH | aimprane | CONTROL | per wn We AT2 (Wa 
CONDITION | | RFA 
Su ce RAD /SEC (Wa), OEGREES 
MEO. BOMBER ELEVATOR $s 546 i75 530 75 
| 
| MED BOMBER, | } 
AIRSPEED = | | ELEVATOR 
isorT/sec | | 
| 300,000L8 
| AIRPLANE ELEVATOR 6 225 429 30 296 
ALTITUDE= | | 
| 
SEA LEVEL | TRANSPORT ELEVATOR 6 23.0 62 24 6¢ 
TRANSPORT RUDDER | 6 566 50 ae9 2 
| MED BOMBER | ELEVATOR 196.7 | 6 
AIRSPEED= | | 
| MED BOMBER, 
| 300,000 LB | | 
ALTITUDE = | AIRPLANE ELEVATOR 5 | 664 318 | 325 216 | 
20,000 FT | | | 
TRANSPORT | ELEVATOR | 6 652 a7 | 216 50 
| | | 
| TRANSPORT | RUDDER | 6 2125 | 36 | 472 os 
| 
— (Crs Cag, + (Chg, 
= —— - 
Cr 
— Gis 


2V —Cig — (Crg)va0 — + 


The spring-tab control is seen to be represented by an inverse 
second-order transfer function, with static gain K, undamped 
natural frequency w,, and damping ratio ¢ which are functions of 
the aerodynamic, kinematic, and inertia properties of the control. 
The spring-tab control’s natural frequency is a most significant 
parameter. Difficulties in obtaining closed-loop stability can be 
anticipated when this frequency is of the order of the natural 
frequencies of the airplane’s short-period modes of motion, as a 
result of the phase lag contribution of the spring-tab control. 
The possibility of this occurrence has been investigated by com- 
paring the undamped natural frequencies of the spring-tab 
control with the natural frequencies of the corresponding short 
period modes of motion, for several typical installations. 

The results of this study are summarized in Table 1. It may 
be concluded that the natural frequency of spring-tab controls 
will usually be sufficiently above that of the airplane’s short- 
period modes to preclude closed-loop stability difficulties from 
this source. In Table 1, the minimum ratio of spring-tab control 
to natural frequency is above 20. For the 300,000-Ib. airplane 
discussed by Phillips, however, the phase lag contribution at the 
estimated closed-loop cutoff frequency (twice the natural fre- 
quency of the airplane’s short-period mode of motion) was almost 
30 degrees, in the low-speed flight condition. This effect is due 
to a high value of spring-tab control damping ratio. This high 
damping ratio was caused by a low effective spring constant 
resulting from low values for C,, and the tab follow-up ratio 6;,. 


The small phase lag contributions at the assumed closed-loop 
cutoff frequencies for the other examples would be unimportant 
in practical installations. 
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ANALYSIS 


MULTICELL torque tube is shown in cross section in Fig. 1. 
The number of cells is arbitrary, but it is assumed that all 
cells are identical. 

Let 6 be the usual torsion stress function, whose directional 
derivative gives the shear stress in the direction of the normal, 
G is the shear modulus, @ is the angle of twist per unit length of 
the tube. St. Venant torsion theory leads to the compatibility 
equation 


v% = 2G0 (1) 
and the boundary conditions 
6 = C; on B,; (2) 


where C; are constants and B, are the boundary curves. Also, 


the torque, 7, is obtained from the relation 


raf 6 dA (3) 
A 


where A is the region enclosed by the outer boundary, Bo. 

The membrane analogy states that 6 may be considered to be 
the deflection of a membrane loaded by a pressure p having a 
surface tension o (force per unit length), where p/o = 2G9. 
The problem of determining 6 is then the problem of finding the 
equilibrium position of a membrane loaded with constant pres- 
sure and having the specified boundary values of Eq. (2). 

An assumed membrane shape near one end of the tube is shown 
in Fig. 2, where we have assumed 6 = 0 on the outer boundary. 

We begin by numbering the cells from one end of the tube, 
so that the membrane deflection over the xth tube is 6z, etc. 

Consideration of the equilibrium of the flat portion of the 
membrane over the (x + 1)st cell leads to the following equation: 


plh = +1 — 82)/telh — 
— 41)/telh + +1)/t/l} (4) 
Replacing p/o by 2G6 and collecting terms, we obtain the fol- 
lowing finite-difference equation: 


6742 2a6, 1 + = —2G6t,l (5) 


* This paper resulted from research carried out under an Engineering 
Study Authorization granted by North American Aviation, Inc., Columbus 
Division. 
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where 
a=1+ (/t-/ht;) 

Linear finite-difference equations with constant coefficients 
yield to substitutions of the form 6; = e in close analogy with 
linear differential equations with constant coefficients. First, 
one must solve the homogeneous equation. Substituting 6, = 
e*, 6241 = eet), etc. into Eq. (5) with the right-hand side set 
equal to zero yields 

— 2ae* + 1] = 0 


Replacing e* by 8, we find the characteristic equation: 


6B? — 208 +1=0 (6) 
This has the roots 
Vat —1 (7) 
The complementary solution of Eq. (5) is, accordingly, 
62 = + (8) 
Assuming a particular solution of the form 6: = K, we easily find 
K = Goht; (9) 
The general solution of Eq. (5) is then 
6: = + + GOht; (10) 


The boundary conditions are that the stress function be zero 
on the outer boundary—i.e., for x = O and x = N + 1, if there 
are N cells in the tube. 

Applying these conditions, we find two equations, Eqs. (11), 
for the determination of C,; and C 

— Goht; (11) 
+ = — GOht; 


Substituting the solutions of Eqs. (11) into Eq. (10), we find 
the equations for the stress function: 


(i = (1 — J 
(x = 1,2,..N) (12) 


The shear stress in the facing sheets is 6,/t;, and this has its 
maximum value at the center of the tube. 

The torque applied to the tube can be calculated from Eq. (3). 
In this case the integration reduces to a summation with 6,4A as 
the summand. AA is the cell area and is equal to h/. Accord- 
ingly, the torque is 


N 
T = >> x 


x=1 


i= +1 +1) Bi? + (aN = +1) Bo (13) 


But note that the sum of a geometric series is given by the formula 


N 
> B7 = B(1 — B)/1 — Bi (14) 


x=1 
Making use of Eq. (14), Eq. (13) can be written in a somewhat 
neater form: 
T = 2G9h*lt; X 


(1 — — BN +1) 


For a one-cell tube (N = 1), Eq. (15) reduces to the familiar 
form given in 


B2(1 — — (15) 


T = 2G6l*h?/[(h/te) + (l/ts)] (16) 


The other extreme, N = ©, is also interesting to obtain from 
Eq. (15). Inspection of Eq. (7) shows that 6; > 1 and & < 1, 
so that high powers of 62 may be neglected, and high powers of 
8, are large numbers. Then it is easy to show that for large N, 
Eq. (15) becomes 


T = 2G9h*lt,{N + [1/(1 — — [62/(1 — (17) 
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If /, the length of each cell, becomes very small while N/ be- 
comes the width of the tube, L, this formula can be well approxi- 
mated by 


T = 2G6h?Lt; (18) 


APPLICATION 


A rectangular four-spar wing of aspect ratio four has been 
tested by Bell Aircraft Corporation and experimentally obtained 
influence coefficients are available.!. From reference 1, we find 
for this wing: 


G = 3.75 X 10° psi, ft. = 0.040 in. 
h = 1.0683 in., ty = 0.063 in. 
= 15/4 in., N=4 
Using these values, we obtain 3;,2 = 6.26, 0.16, from which 


Eq. (15) yields a torsional stiffness of 7.31 X 10° in.-lb. per 
radian. Applying unit loads at the leading- and trailing-edge 
spars at the wing tip (to form a couple) and then averaging 
wing rib rotation angle as obtained from the experimental in- 
fluence coefficients for these loads, we find an experimental 
torsional stiffness of 7.30 X 10° in.-lb. per radian. This excellent 
agreement seems to be adequate confirmation of the utility of 
classical torsion analysis to even very-low-aspect-ratio wing 
structures. 
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SUMMARY 
on empirical formulas are given for pressure coefficient 
and pressure-drag coefficient for axisymmetric bodies with 
hemispherical noses in both subsonic and supersonic airflow, at 
zero angle of attack. 
Comparison with existing test results shows excellent agree- 
ment for all speeds. 


INTRODUCTION 

It is known that both the pressure coefficient of a hemisphere 
and the pressure-drag coefficient derived therefrom, either by 
using the Newtonian theory or by using its modified version, do 
not agree very well with those obtained from test data, apart 
from the fact that these approximations cover only the super- 
sonic and hypersonic flow regimes. 

Other approximations for pressure drag at all speeds have 
been suggested in references 1, 2, and 3 with the main objective 
of facilitating calculation of trajectories with computing ma- 
chines. Recently, the present writer obtained two rather simple 
empirical expressions, one for pressure coefficient, the other for 
pressure-drag coefficient. They correlate the test data very 
well in both the subsonic and supersonic range. 

These fits are simple to apply, hence can be used as design cri- 
teria and test guide. Since they all contain an arc-tangent 
term, they will be called arc-tangent fits. 
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TABLE 1 
Pressure Coefficient Calculated by Different Methods 
M 0 1/2 3/4 1 3/2 2 3 © 
Newtonian 
Modified -- -- -- 1.83 | 1.83 1.83 | 1.83 
M= 
Newtonian 
Modified 
4 |M#a@ -- -- -- 1.00] 1.46 | 1.62 /1.74 | 1.83 
First Arc 
Tangent 1.00) 1.10 |1.21 1.35] 1.60 | 1.70 1,85 
Fit 
1.00] 1.10 1.2 
1.30] -- 1,66 1,75 1.83 
2.0 
16 
- —-— Modified Newtonian for M= © 
Modified Newtonian for Mé¢ @ 
—— First Arc-Tangent Fit 
12 
8 ° / © |NOL unpublished (see Reference 9 ) 
| Reference 5 
0.8 © | Reference 6 
o> | Reference 7 
OC | Reference 8 
04 
% 3 4 6 7 
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Fic. 1. Maximum pressure coefficient for a hemisphere-cylinder. 


TABLE 2 
Comparison of Pressure Drag Coefficient, Calculated and 
Experimental 
M o | 1/2 | 3/4 1 | 3/2 2 3 4 co 
Newtonian 
Modified -- | -- | -- | -- |.913 }.913 | .913] .913 | .913 
M= 
Newtonian 
Modified -- | -- | -- | -- |.730 |.800 | .807] .900] -- 
M# 
ra) 
° | Second 
Arc .12 | .20 | .26 | .36 |.600 |.750 | .850].890 | . 960 
Tangent 
Fit 
Test 20 22 
-23 | .34 |.650 |.750 | .830].875 | -- 


a | 4 


—-— Modified Newtonian for M=00 
Modified Newtonion for 
—— Second Arc-Tangent Fit 


j 
C ——— | 
© | NOL Unpublished (Reference 9) 
| _Reference 5 _| 
od Reference |i 
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MACH NUMBER 
Fic. 2. Pressure-drag coefficient for a hemisphere. 
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(1) Pressure Coefficient as a Function of Mach Number 

Denoting the free-stream Mach number by / and the maxi- 
mum pressure coefficient by Cpmar, the modified Newtonian for- 
mula for Cpmaz given by Oliver‘ reads: 


Cpmaz = (11/6)[1 — (5/11)-(1/.M?)] (1) 


where the ratio of specific heats for air istaken as 1.4. This equa- 
tion is an approximate expression for Cpmaz as a function of J/. 
A better approximation is provided by the first arc-tangent fit. 


Comaz = 1.85 — (1/7)are tan [1.96 (1 — M)] (2) 


To show this, values of the pressure coefficient are calculated 
by using (1) modified Newtonian theory for M = o, (2) modi- 
fied Newtonian theory for M # ~o, (3) the first arc-tangent fit. 
Values obtained from these methods are compared with the test 
results from references 5-8, and the NOL data given in reference 
9. The results are shown in Table 1. 

The curve in Fig. 1 at the end of this report shows clearly that 
the arc-tangent fit can be used as a design criterion for both 
supersonic and subsonic range. 


(2) Pressure-Drag Coefficient as a Function of Mach Number 


Ordinarily the pressure-drag coefficient is either obtained by 
integrating the pressure coefficient given by the modified New- 
tonian theory or calculated directly from test data by integration. 
Here we will avoid the complicated integration given in standard 
textbooks and use the second arc-tangent fit instead: 


Cp = 0.46 — (1/m)arc tan [1.61(1.2 — M)] (3) 


where M denotes the Mach number and Cp the pressure-drag 
coefficient. 

Table 2 shows good agreement between the second arc-tangent 
fit and test results as compared against those computed from 
Newtonian theory and its modified version. 

Table 2 shows that the values of pressure-drag coefficient cal- 
culated by the second arc-tangent fit are in good agreement with 
test results in the transonic and supersonic range. In the sub- 
sonic range the maximum deviation from test data is 16 per cent 
for M > 1/2. This means even in the subsonic region the second 
are-tangent fit closely predicts the pressure-drag coefficient. 
The curve shown in Fig. 2 again indicates that the simple arc- 
tangent fit agrees quite well with test data in both subsonic 
and supersonic regimes. 


CONCLUSION 
Both the pressure coefficient and the pressure-drag coefficient 
for an axisymmetric body with a hemispherical nose predicted 
by the arc-tangent fits are in very good agreement with existing 
test results, while the Newtonian theory and its modified ver- 
sion, both of which are limited to the supersonic range, give re- 
sults which are not in good agreement with experiments in that 
region. 
In addition to their simplicity, the arc-tangent fits are unique 
for all speeds. They not only apply in the supersonic range 
but are also valid for subsonic speeds. 
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— PROBLEM of minimum time-to-climb of an airplane 
which was analyzed by the author in 1957! has recently been 
reconsidered by Theodorsen.? Without directly including the 
equation for the weight in his treatment, Theodorsen has brought 
considerable simplification to the numerical analysis of the 
problem. An alternative method for solution (which was im- 
plied in the author’s formulas’) is given below. 

If we take the height / as the independent variable, the time- 
to-climb of an airplane in the (x, h) plane is given by 


tn = sin (1) 


where V is the speed, and y the positive inclination of the path 
with the horizon x. The problem is to find the path of an air- 
plane of given characteristics which corresponds to the shortest 
time between two prescribed points such that V and y at these 
points have given values. This problem has been analyzed 
broadly in the work of the author. Let us briefly recall the re- 
sults. If the relative decrease in the weight of the airplane dur- 
ing the climb is neglected, and if the distance x to be covered is 
not specified, then the equations governing the problem are? 


V’ = N(Cr — Cop) — (g/V) (2) 
(N/V)CxL — (g/V?2) cot y 


and 


Mi’ + + Aare + As 
do’ + + Bede + Bs 
— (1/V)(2/M) = 0 
where Cr, Cp, and C, are conventional coefficients for thrust 
T, drag D, and lift L, 
Cr T/[(1/2)p V2S], ete. 
N = gpSV/2W sin y 
A, = (N/V)(Cr — Cp) + N(0/0V)(Cr — Cyn) + (g/V?) 


0 
0 (3) 


A, = (22/V*) cot y, A; = (1/V? sin y) 
B, = —(Cr — Cp)N cot y, (4) 
B, = —(1/V)[NCz cot y — (g/V sin? y)] 


= 


and dy, 2 functions of altitude, are the Lagrange multipliers 
corresponding to Eqs. (2); the prime indicates derivatives with 
respect to h. One may, of course, write 0/0V = (0/0M) 
(1/a), where a = speed of sound. Eqs. (3) may be compared 
with Eqs. (16) of Theodorsen. 

The d’s can now be eliminated from Eqs. (3), and a single equa- 
tion for £, the actual unknown of the problem, results. This equa- 
tion can be obtained as follows. From the last equation in Eqs. 
(3) we immediately obtain 


B; = (cot y/V sin y), 


do’ = EVA,’ (5) 
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The first two equations of Eqs. (3) can now be rewritten 


+ (A; + + As = 0 


EVA! + [(EV)’ + Bi + + By = 0 (6) 


Hence, 


— + A(EV)’ + — A, Bs] + 
(A, — + AXEV)? — (EV)’ — B] = 
(d/dh){(B; — AsxEV)/{(A, — BEV + 

A(EV)? — (EV)’ — Bi} (7) 


which is a differential equation for {V = X of second order and 
third degree. This may be written out as 


(Bs — AsX)X" + GX’? + GQX'X + 
+ GX? + OX? + +G =0 (8) 

For the sake of brevity, we shall not give here the values of co- 
efficients C,, C2, . . . Cz, which are simple functions of A’s, B’s, 
and their first derivatives with respect to h, and which are easy 
to evaluate from Eq. (7). The integration of Eqs. (2) and (8) 
involves four arbitrary constants which can be determined by our 
four boundary conditions, namely: V = Vo, y = yo ath = 0; 
V = Va, y = yu ath = H, where Vo, yo, Va, yu, and H are 
given. Under these conditions, the problem has a unique solu- 
tion to be determined by a step-by-step numerical integration. 

We have tacitly assumed that the coefficient of X¥” does not 
The vanishing of this coefficient occurs when 


= cot (9) 


The poiats which satisfy Eq. (9) constitute singular points of 
solution and should be avoided in flight, if they exist (the coeffi- 
cient of X” changes sign in their neighborhood). Such a singular 
point appears, for instance, at the end of the path of an unre- 
stricted climb when there are no assigned values for Vy and yx 
(see reference 1). 

The contributions of both Theodorsen and the author thus con- 
clude a problem of great military significance. 


vanish. 
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SYMBOLS 
a = sound velocity, isentropic 
e = specific internal energy 
h = specific enthalpy 
Mach number = 2/a 
~p = pressure, ~/Po = pressure in atmospheres 
v = velocity component normal to shock front 
w = velocity, defined in text 
a = incident angle 
y = isentropic exponent = pa*/p 
p = density 
6 = flow displacement angle 
Subscripts 
1 = flow upstream of shock 
2 = flow downstream of shock 


6 pe MEASUREMENT Of Mach angle is a useful means of deter- 
mining some of the flow characteristics on a body in hyper- 
sonic flow. In particular, it has proved valuable in the deter- 
mination of pressure distribution in shock-tube experiments.! 
The relation between Mach angle and Mach number when the 
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gas (air) is in dissociation equilibrium may be derived in a simple 
form comparable to the usual perfect-gas Mach-angle equation. 

From the conservation of mass, momentum, and energy across 
a normal shock, we may obtain the Hugoniot relation in the 
following form :? 


— hy = [(pe p1)/2p1] (pi/p2)] (1) 
from which 

b2/p2 = (pr/p2){h2 + [he — (p2/p2)] — — 

— + (p1/p2)} 
Or, since h = e + (p/p), 
p2/p2 = (pi/p2)[(he + e2) — (ty + + (P1/p2)] (2) 
Now, define 
Ki = (pi/p1)/(y + Ka = (p2/p2)/(he + (3) 
Again, from the normal-shock relation, 
= 1 + yM?[1 — (4) 
where y is defined as: 
y = pia,?/p, = (Oln p/d In p)s 
Inserting Eq. (3) and Eq. (4) into Eq. (2), we obtain the quadratic 
relation for p;/ps: 
(ps/p2)* — [((K2 + 1)(1 + (01/2) + 
K2[1 + (1 + Ki™')y7M~*] = 0 
The solution of this equation, after considerable algebraic 
manipulation, may be written as: 
pi/p2 = (1/2)(Ke + 1)(1 + + 
(1/2){ — Ke) — (1 + + 
Note that if the gas is perfect, then K2 = Kj; if positive root is 
considered, then p;/p2 = 1, which is the trivial solution, hence 
negative root is chosen. 

Now, as in the derivation of the Mach angle for a perfect gas, 
consider an oblique shock—i.e., impose a velocity perpendicular 
toa normal shock. From Fig. 1 and continuity, we have 

tan (a — 0)/tan a = = pi/pe 
Now let 6 — 0, then Kz — K, and, from Eq. (5), we have, upon 
substituting G = M, sin a for M, 
(i~ (©) 


or 


= (1 + K)/[y(1 — K)) 


Again, if we ccnsider a perfect gas, then K = (y — 1)/(y + 1) 
and 


(1+ K)/[y1 K)] =1 


Fie. 1. 


4 COMPUTED FROM REF 3 + 
105 COMPUTED FROM REF 4 — 
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Fic. 2. 


Hence, we take the positive root of Eq. (6) or 
sin a = (1/M,)[(1 + K)/7(1 — K)]'/72 (7) 


The factor [(1 + K)/7(1 — K)]'/2 has been computed approxi- 
mately using the real-gas curves of references 3 and 4 and is 
plotted in Fig. 2 as a function of pressure with temperature as a 
parameter. 
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N THE WHIRLING unbalanced shaft there are three points of 

importance: the center of the bearing, O, the center of per- 
cussion of the shaft, W, and the center of gravity of the disc, S. 
Neglecting friction and gravity and without oscillations, these 
three points lie on a straight line. The straight line OWS (be- 
low critical speed) or OSW (beyond the critical speed) in that 
order is supposed to describe a stationary motion about the 
point O. 

But the surfaces of the disc are invariably subject to friction 
from the environment in which they operate. The effect of 
environmental friction can be separated into a frictional moment 
and a single force. The frictional moment is overcome by the 
useful tangential force (or by an external turning moment trans- 
mitted through the shaft) and need not be considered further. 
The single force acting at the eccentric center of gravity on the 
other hand disturbs the equilibrium, or the stationary motion of 
of the three points, O, W, and S only in the neighborhood of the 
critical state. In other words, environmental friction destroys 
the collinearity of OWS as the shaft is being speeded through the 
critical state. 

Assuming that the disturbance impressed by environmental 
friction force on the stationary motion of the shaft (OWS or 
OSW) is a small perturbation, this note proposes to find out by 


*The author’s thanks are due to Professor A. Ramachandran for his interest 
in this project and his constant encouragement. 
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the method of small vibrations whether the disturbance increases 
or decreases and, if so, under what conditions. 

It is assumed that the frictional force is opposed to the ab- 
solute velocity of the center of gravity and has a magnitude 


R= (1) 


where uw is a (not easily estimable) coefficient of environmental 
damping and m, the mass of the disc, is introduced for dimensional 
purposes. This assumption is well within the experience of rota- 
tion in a viscous medium. 

The motion of the disc is considered relative to a rotating co- 
ordinate system (Y, Z) rotating with an angular velocity a, 
> about the point O. 

With @ as the relative angle of rotation, the equations of 
motion of the ‘‘disturbed motion”’ are (see Fig. 1) 


mi = myw? + 2muyz — aly — e,) — pmyar? | 


mz = — — — e,) — (2) 
I6 = —al(y —e,)e, + — e,)e, 


Treating y, z, and @ as quantities of small magnitude and neg- 
lecting the products and powers of small magnitude, Eqs. (2) 
reduce to 


— By — 2mayz = ew,? (3a) 
2 — Bs + — ew,20 = O} (3b) 
e6 — (e2/k?) w,22 = 0 (3c) 
where 
B = — — w,?] 
and 


I = mk? 

Since the eccentricity e at the actual disc used in practice always 
remains very small compared to its radius of gyration k, say 
e/k = 107%, then e?/k? is the order of 10-6 and can be neglected. 
Eq. (3c) gives for 0 


which disappears identically as soon as it is assumed that wy 
is the mean velocity of rotation of the disc and if the Y-axis is 
located in such a way that at ¢ = 0 it lies in the direction WS = 
é. 
By a transformation y = » + 6 and by choosing 6 suitably, 
Eq. (3) reduces to 
7 = Bn — 2mayz = 0 (4) 
— Bs + 2mann = O } 
The trial solution 
n = ae“, = be 
leads to the frequency equation 
Af + — B] + B? = (5) 
It can be shown that for small values of the difference (a, — 
w;), stability is attained as soon as 
(wo? — > (6) 
Neglecting environmental friction but taking into account 


e*/k?, Stodola,' from entirely different considerations, established 
the condition for stability beyond the critical speed as 


[(wo? — > 4(€/k?) (7) 


Eqs. (6) and (7) are independent criteria of stability. The 
former takes into account the friction of the environment and 
neglects e?/k?, while the latter neglects friction but takes into 
account e?/k?. 

To decide which one of the above two inequalities is satisfied 
earlier, it is necessary to find out whether 


(8) 


To decide this, the deflection of the shaft at the critical speed 
has got to be determined. From Fig. (2), the deflection 7,,,, for 
small y can be shown to be 
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Fic. 1. Forces and apparent forces in the rotating system. 


FIs. 2 
Oo = bearing center 
Ww = shaft center (at the disc) 
8 = gravity center 
Ws =e = eccentricity of the shaft; ¢y, ¢2 components of (y, 2) di- 
rection 


7 = displacement of the shaft 

BY = displacement of the center of gravity of the disc 

F = outward centrifugal force = mw*y 

R = frictional force opposed to the absolute velocity of the center of 
gravity of the disc = umy w* 


P = elastic force of the shaft = ma ,?*r 
wk = critical speed of the shaft if no friction is present 
wo = angular velocity of the disc in steady rotation 


tan y = R/F =u 
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?mar = e/p (9) 


sin y =tany =u (10) 


ESTIMATION OF pu 
For a disc whose e and k are known, assume as a first approxi- 
mation of u from Eq. (8) is given by the equality 


p= W4(e2/k?) (11) 


With this value of 4, compute the maximum deflection 7,7 from 
Eq. (9). Compare this value of %maz with that obtained by a 
simple experiment on the disc whose e and k are as given above. 
If experimental value agrees with the computed value, then the 
assumption of Eq. (11) hoids. If the computed 7mgr is greater 
than the observed 7mqz, then w as given by Eq. (11) is an under- 
estimate and, hence, one may conclude that yu is positively greater 
than that given in Eq. (11). Furthermore, y» according to 
Eq. (10) must be less than 1. Therefore, » should lie between 


4(02/k?) and 1, that is, 
V 4(e2/k?) <1 (12) 


Plotting 7mar Vs. u from Eq. (9), the curve as shown in Fig. (3) 
is obtained and the value of yp lies in the shaded region. 

By computing the area of the shaded portion of the curve be- 
tween the limits V 4(e2/k?) and 1 and dividing it by the base, a 
mean value of 7mgz is obtained. A mean value of yu is obtained 
by back substituting this value of 77 in Eq. (9). 

For example, in a test reported by Stodola? on a rotating disc 
with eccentricity e = .OL” and radius of gyration k = 5”, the 
maximum deflection at the critical speed w) = w,; found experi- 
mentally was 0.0375”. 

Applying the procedure outlined above, the value of uw for this 
case was found to be 0.27. 


REFERENCES 
1 Stodola, A., Steam and Gas Turbines, Vol. II, p. 1133, Peter Smith, 
New York, 1945. 
2 Ibid, Vol. 1, p. 429. 


A Rapid Method for Estimating the Shear 
Stress and the Separation Point in Laminar 
Incompressible Boundary-Layer Flows* 
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Department of Chemical Engineering, University of California, 
Berkeley, Calif. 

November 13, 1959 


SYMBOLS 
cy = dimensionless shear coefficient at the wall = 70/pl’ ~? 
L = characteristic length of the surface 
Ue = characteristic velocity of the system 
U = velocity at the edge of the laminar boundary layer 
U1 = U/JUe@ 
u = x-component of the velocity 
um = 
v = y-component of the velocity 
= Re 
x = distance along the surface from the leading edge 
x1 = x/L 
y = distance normal to the surface 
(y/L)V Re 
v = kinematic viscosity 
p = density 
70 = shear stress at the wall 
Re = Ua@L/v, the Reynolds number 


_ PURPOSE of this note is to present a rapid method for 
estimating the shear stress and the separation point in laminar 


* This work was supported by a grant from the National Science Founda- 
tion. 
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TABLE 1 
(The values for F,,”"(0) have been obtained from reference 6.) 


n F,,"(0) [(n + 2)/4]2/3F,,"(0)4/3 
—0.181 0 0 
—0.178 0.0552 0.0124 
—0.174 0.0857 0.0224 
—0.165 0.129 0.0388 
—0.148 0.191 0.0658 
—0.131 0.240 0.0898 
—0.095 0.319 0.133 
—0.049 0.400 0. 163 

0 0.470 0.230 

0.051 0.531 0.275 
0.105 0.587 0.320 
0.222 0.687 0.409 
0.353 0.775 0.500 
0.500 0.854 0.592 
0.667 0.928 0.691 
0.857 0.996 0.795 
1.000 1.039 0.869 
1.333 1.120 1.0380 
2.000 1.233 1.322 
3.000 1.336 1.707 
>3 1.2V2n/(n + 2) 0.80 2/3 


incompressible boundary-layer flows past arbitrary surfaces. 
We begin with the basic equations in their well-known dimension- 
less form: 

)+ )= ( 1/2 (dU 2/dx; ) + (1) 


and 


+ = O (2) 
Let us suppose, for the moment, that 
U;? = a,x," + (3) 


where a, m, and are arbitrary constants, with aj, a2 > 0. 
In view of Eq. (3), it is apparent that c;, the dimensionless shear 
coefficient at the surface, is of the form 


V Re cy = =0 = Z(X1, a1, a, M1, (4) 


where Z is a function of the indicated parameters. 

Now, it is obvious that, although the exact evaluation of Z 
is in general a matter of considerable difficulty, its asymptotic 
form, as either a; > © or ay > ©, is already known. Thus, 


Z( an, a, > + 2)/4] "(0) (5) 


as a,j > ©, with, of course, a similar result when ax —~ ©. The 
parameter F,,"(0) can readily be obtained from Hartree’s numeri- 
eal solution? ® of the Falkner-Skan equation® and is shown in 
Table 1 for various values of n. 

It is reasonable to suppose next that, for intermediate values of 
a and a, Z may be estimated by a simple interpolation formula. 
We choose the simplest possible expression 


ym + ( had m (6) 


where m is a constant. 

However, since Z must reduce to the exact solution, Eq. (5) 
with a = (a; + a2), when = m, it is clear that m must be set 
equal to 4/3. 

In general, by suitably choosing the parameters a; and nj, the 
latter not necessarily integers, we may represent U’,? by 


k 
i=1 
to any desired degree of accuracy. Then, in view of Eqs. (5) and 
(6), we are led to believe that the dimensionless shear coefficient 
at the surface may be computed approximately from 


k 3/4 
V Re cy = {> + 2)/4)? Bain” 


i=1 
(8) 


Strictly speaking, this result is rather empirical, especially since 
we shall proceed to make use of it even if some of the aj’s are 
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TABLE 2 
xX) 
(separation) (separation) 
Eq. (9)| exact reference 

(1 — x,)? 0.149 0.126 (2) 
(1 — x,?)? 0.306 0.290 (3) 
(1 — x,‘)? ().485 (3) 
(1 + x)" 0.189 0.160 (2 
(x; — 6.29 10> 4x,3 — 

4.62 X 10~5x,5)? 6.70 6.87 (2) 
(sin x, )? 1.93 1.93 (5) 


negative. It appears, however, to give results which are in good 
agreement with those obtained from exact solutions. Further- 
more, it is quite probable that this general procedure could also be 
used successfully for estimating both the displacement and the 
momentum thickness which are usually of interest in boundary- 
laver calculations. 

In order to test the accuracy of Eq. (8) at the most sensitive 
spot, the separation point has been calculated for a number of 
simple potential-flow distributions, which belong to a special 
class of functions Uj(x;) for which exact series solutions have 
been made available.2-> The comparison is made in Table 
2, from which it is apparent that in every case the 
empirical interpolation formula gave answers in reasonable 
agreement with the exact results. It was also found that, 
whether the pressure gradient outside the boundary layer was 
favorable or not, the shear stress could be computed from Eq. (8) 
with an error which for all practical purposes could be ignored. 

It is believed, therefore, that, for estimating the shear stress 
and the separation point in those numerous cases where exact 
solutions have not as yet been developed, Eq. (8) might be used in 
preference to the usual approximate methods of boundary-layer 
theory,': °~7 which are not only far more complicated than the 
simple procedure just described but, unfortunately, no more accu- 


rate or reliable. 
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2 Gortler, H., A New Series for the Calculation of Steady Laminar Boundary 
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3 Gortler, H., and Witting, H., Zu den Tanischen Grensschichten, Ostr. Ing. 
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mate Treatment of the Equations of the Boundary Layer, Proc. Cambr. Phil 
Soc., Vol. 33, p. 223, 1937. 

Schlichting, H., Boundary Layer Theory. McGraw-Hill Book Co., Inc., 
New York, 1955. 

€ Smith, A. M. O., Rapid Laminar Boundary-Layer Calculations by Piece- 
wise Application of Similar Solutions, Journal of the Aeronautical Sciences, 
Vol. 23, No. 10, pp. 901-912, October, 1956. 

* Wieghardt, K., Uber einen Energiesatz sur Berechnungen Laminarer 
Grenzschichten, Ing. Archiv., Vol. 16, p. 231, 1948. 


Computer Solution Techniques for 
Two-Dimensional Supersonic Flow Problems; 


William P. Harrison, Jr.* 
Experiment Incorporated, Richmond, Va. 
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SYMBOLS 
Mo = free-stream Mach number 
Momin= free-stream Mach number, below which leading-edge shock wave 
will be detached for given value of 7 
M, = local Mach number at point P 


+ The research program under which this work was done was supported 
by the Bureau of Ordnance, Department of the Navy, under Contract 


NOrd-18030. 
* Presently on academic leave of absence at Virginia Polytechnic Insti- 


tute, Blacksburg, Va. 
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K = 1) 

a = angle of attack 

7 = ratio of specific heats for air 

é = slope of surface at leading edge 

by = slope of surface at point P 

7] = flow turning angle at leading edge 

nmax = flow turning angle above which leading-edge shock wave will be 
detached for given value of VW» 

6, = wave angle of leading-edge shock wave 


INTRODUCTION 
princes involving two-dimensional supersonic flow around 
concave corners and over convex surfaces are often solved 
on digital electronic computers, and the computational methods 
employed largely determine the accuracy of final results as well 
as total machine time required. 


METHOD DEVELOPMENT 
In Fig. 1, free-stream flow is shown approaching convex sur- 
face OA at an angle a, such that the flow first must negotiate 
the concave corner (6; — a@) before expanding over surface OA. 
For accurate machine calculation of the leading-edge flow condi- 
tions it is convenient to use the wave-angle cubic equation of 
reference 1, assuming a, Wy, and 6; are known quantities. With 


substitutions 7 = 6; — a and 8 = sin? @,, this cubic can be 
written 

bB2?+ (la) 
where 
b = —[(y sin? n + 1)Me? + 
= | [0.25(y + + y — 1) Mi? sin? + (1b) 


2M,? + 14/4 
d = —(cos? n)/ 
Methods involving iterative procedures are popular for the 
solution of Eq. (la). The exact solution, however, is naturally 
more direct, more accurate, and less time consuming. With 
regard to the roots, 8), 82, and §;, there exist three possible solution 


combinations: 


(a) One real root and two conjugate imaginary roots. 
(b) Three real, unequal roots. 
(c) Three real roots of which at least two are equal. 


Solution (a) occurs when My < M.,i,@ n or, to express it dif- 
ferently, when 7 > 9mar@ Mo. For this detached-shock condi- 
tion, the real root which can be extracted from Eq. (la) has no 
physical significance. 

Solution (b) is the root combination most often encountered. 
As noted in reference 2, the smallest root corresponds to an en- 
tropy decrease and should be disregarded. If the magnitudes 
of the three roots are such that 6, < 62 < 8;, then only 8: and 8; 
remain for consideration. The largest root, 83, gives the strong 
shock-wave angle; thus, middle root 8. is desired. The equation 
for B. can be expressed 


{ [cos (6/3) — V3 sin / 


where 


SHOCK~__ Z-LEADING - EDGE 
WAVE TANGENT 
J 8p 
= 
P(x,y) 
Ow 


Fic. 1. Expansion surface with sharp leading edge in super- 
sonic flow. 
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Fic. 2. General shape of the variation of the function g( 17) 
with Mach number. 


H, = —(b?/3 — c) } 
H, = —[(b/3)? + H,](6/3) + d > (2b) 
@ = V + 1] 
with the requirement that ¢ always be taken as principal value 
of the arc-tangent fuaction. 

The third root possibility, mentioned under (c), is a special 
case of (b) and occurs only for that particular combination of 
M, and n which results in a limiting shock-attachment condition 
at the leading edge of surface OA. That is, Wy = Monin@ 
n OF 7 = Nmar@ My. The two equal values of @ which evolve are 
the wave angle for both strong and weak shocks, since the two 
waves coincide at this condition. The only computational con- 
sequence of this is that @ goes to zero; thus, Eq. (2a) still holds. 
Again, the one remaining root, which will be smaller than the 
two equal roots, is physically insignificant. 

With 8 determined, the local Mach number behind the leading- 
edge shock wave can be calculated using Eq. (132) of reference 2: 


— (y — — + 2] 

At any local point, P, on expansion surface OA, the familiar 

Prandtl-Meyer equation 

y = M? — 1/K) — VM? (4) 
and its first derivative 

dv/dM = 2V M? — 1/M((y — 1)M? + 2] (5) 

can be used to determine local Mach number, from which other 
flow quantities can be computed. The development to follow 
presents a simple machine method for the rapid calculation of 
M,, the local Mach number at point P. 

In the most usual case, profile curve OA is fully defined so 
that the local slope, 6,, is known at point P. Substituting the 
functional notation 

vy = f(M) (6) 
and 
dv/dM = g(M) (7) 
it can be stated that 
— = f(M,) —f(M) = 


for either an upper- or lower-surface profile, and 


ls, ~ & 


Mp 
— 61! e(M)dM (8) 
M 


M 


The general shape of the function g(J/), with a value of 1.4 
for y, is shown in Fig. 2. Considering the distribution of area 
under this curve and Eq. (8), it becomes apparent that for a first 
approximation 


Mn = M; 
pl 1 + 2(M;) 


(9) 
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For a second approximation, then, 
Mp: = Mn + {LF — (10a) 
where 
F = + — (10b) 


This procedure can be continued to give as close an approxi- 
mation to J, as desired, since the process is always convergent. 
Expressed in more general terms, 


M, = } (11) 
Mon = Mon-1 + {(F — f( Mon—-1)] /2( Mopn-1)} } 
(42) 


where F is given by Eq. (10b) and 1/,, is given by Eq. (9). 


REFERENCES 
1 Thompson, M. J., A Note on the Calculation of Oblique Shock-Wave 
Characteritsics, Journal of the Aeronautical Sciences, Vol. 17, No. 11, p. 766, 
November, 1950. 
2 Ames Research Staff, Equations, Tables, and Charts for Compressible 
Flow, NACA Rpt. 1135, 1953. 
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Ends 
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I A RECENT NOTE,! Neffson outlines a procedure for computing 
the end moments and related quantities for a beam-column 
span of constant section with known (or assumed) elastic rota- 
tional restraints at its ends. The proposed procedure is need- 
lessly complex and tedious because no advantage was taken of 
beam-column formulas and tabulated functions that have long 
been known in the aircraft industry and are most extensively 
piesented in reference 2. It is true that the case discussed by 
Neffson is not covered in Airplane Structures, and the purpose 
of this note is to indicate how the formulas and tabulated func- 
tions in that volume can be utilized to handle Neffson’s problem. 
In Art. 14:3 of reference 2, a three-moment equation is de- 
veloped for beam-columns by equating expressions for the slopes 
(or rotations) of spans meeting at an intermediate support. In 
order to improve the parallelism of the result with the tradi- 
tional three-moment equations for continuous beams without 
axial loads, all terms of the basic continuous beam-column for- 
mula are multiplied by 6 EZ. The resulting three-moment equa- 
tion for a continuous beam-column subjected to a uniformly dis- 
tributed transverse load in each span is presented as Eq. 14:11. 
The method of adapting it to other types of transverse load is 
covered in Art. 14:5 of reference 2. 

Consider the beam-column shown in Fig. 1, assuming that the 
reaction moments, M4 and Mz, are produced by elastic restraints 
at supports A and B and are directly proportional to the rotations 
of the member at those supports. In view of the method used 
for developing Eq. 14:11, we can write the two equations 


= —2L8M4 — LaMe + 0.25L3 yw (la) 
6El0p = LaMas + — 0.2513 yw (1b) 


If the beam-column were fixed at A and B, 64 and @g would be 
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zero, and Ma, and Mg would be the fixed-end moments for the 
given loading. If we assume 64 = 6g = OU and solve the above 
equations for 1/4 and Mz, we get 


Ma = Mp = 0.25L?yw/(268 + a) 


This agrees with the fixed-end moments for uniformly distributed 
side load, as given in Art. 14:11, when proper account is taken 
of the sign conventions used. 

The restraint provided by the supports at A and B can be ex- 
pressed quantitatively by the ‘‘flexibility factors” 


Fa = 6EI0@4/Ma 
= 


The plus sign is used for F4, since a counterclockwise (or +) 
value of 64 is recuired to develop a positive value of M4. The 
minus sign is used for Fxg, since a clockwise (or — ) value of 62 is 
required to develop a positive value of Mg. The factor 6EI is 
inserted for computational convenience and because of its use 
in developing Eqs. (la) and (1b). It is assumed that numerical 
values of F4 and Fg can be obtained from the properties of the 
supporting structure and EJ of the beam-column. Both are 
inherently positive, like length and EJ. 
Then, from Eq. (la) and (1b), 

Fs = +6E/04/Ma = —2L6—LaMg/Ma + L'yw/4Ma (2a) 
Fx = —6EI0g/Mg = —LaMa/Mpg — 2L6 + (2b) 


whence 


+ 2L8)M4 + = L*yw/4 (3a) 
+ (Fe + 2L6)Mz = L*yw/4 (3b) 


For any specific beam-column one can compute L/j trom the 
given values of P, L, E, and J and use Table 14:2 to find the 
values of a, 8,and y. For other types of side load, suitable terms 
for the right-hand sides of Eqs. (3a) and (3b) can be obtained 
from the columns of ‘‘Terms for Three-moment Equation”’ in 
Table 14:3, by multiplying the expressions listed in that table 
by J. The term to insert in Eq. (3a) should be taken from the 
“Right Bay,” and that for Eq. (3b) from the ‘‘Left Bay’’ column. 

Thus, for the case of a concentrated side load shown in Fig. 2, 
the equations for computing the reaction moments would be 


(Fa + 228)M4 + LaMp = 
6W7?{ [sin (b/j)/sin (L/j)] — (b/L)} (4a) 


LaM, + (Fe + = 
6W7?{ [sin (a/j)/sin (L/j)] — (a/L)} (4b) 


Similar equations for more complex loadings can be constructed 
by the method of superposition of loading terms on the right- 
hand sides. 

As can be seen from Eqs. (3), if the span is symmetrical with 
respect to both loading and support stiffness, the reaction mo- 
ments will be equal and can be computed from a single equation. 

If a beam-column were of constant section in each span, and 
continuous over three or more supports and elastically restrained 
against rotation at each of them, three-moment equations could 
be developed by suitable manipulations and combinations of 
equations like Eqs. (2), It would be much simpler, however, to 
use moment distribution along the lines presented in Art. 14:11 
of reference 2. 

In this event, each support would be considered as an addi- 
tional member with an absolute stiffness 17/6, and the unbal- 
anced moments at each support would be distributed to the two 
beam-column spans and the support itself in proportion to their 
relative stiffnesses. In practice, it would usually be advisable to 
use relative support stiffnesses of 17/46, since the stiffness co- 
efficients of p. 121 of reference 2, when multiplied by J/L, give 
only one quarter of the ratio of moment applied at one end of a 
beam column to rotation at that end. 

Fixed-end moments, carry-over factors, and stiffness factors 
for the beam-column spans would be computed by the methods 
of Art. 14:11. If only moderate precision were desired, they 
could be obtained by the use of Figs. 14:12-14:22 inclusive. 
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Fic. 2. 


More precise values could be obtained for most, if not all, of the 
quantities of interest by using Eqs. 14:27—14:33 or Tables 14:8 
and 14:9. 

Occasionally, it may be desired to use more precise values of 
fixed end moment for a loading which is not adequately 
covered by those formulas and tables. Then, it would be found 
appropriate to write pairs of Equations like Eqs. (la) and (1b) 
with 64 and 6g set equal to zero and to solve them for the desired 
fixed end moments. 

Since practically any transverse loading likely to be encountered 
in practice can be closely approximated by a combination of the 
loadings for which three-moment equation terms are given in 
Table 14:3 of reference 2, the method outlined here provides a 
simple straightforward method for the analysis of beam-columns 
with supports elastically restrained against rotation. 
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A Note on the Relation of Biot’s Method in 
Heat Conduction to a Least-Squares 
Proceduret+ 


Stephen J. Citron* 
Institute of Flight Structures, Columbia University, N.Y. 
November 18, 1959 


_ PURPOSE of this note is to show that a relation exists be- 
tween Biot’s method in heat conduction!~* and the method of 
least squares. In doing so, a means will be given for judging the 
relative accuracy of several approximate solutions determined by 
Biot’s method. 


Brot’s METHOD 


For a system within the region 2 bounded by the surface 2, 


the method requires that the function 


v= f +e an (1) 


be made a minimum by varying the heat-flux vector (H) while 
holding the thermal force (V7) constant. To complete the for- 
mulation, it is necessary to require that energy be conserved in the 
system—i.e., 

v-H = — pT (2) 
The use of this relation with the function ¥ has been shown by 
Biot to be directly equivalent to satisfying the Fourier heat con- 
duction equation. 


+ The author would like to thank Professor Bruno A. Boley of the Institute 
of Flight Structures, Columbia University, for the advice and assistance 
rece‘ved during the course of this investigation. The work is part of a 
project supported by the National Science Foundation under grant number 
NSF-G5948. 

* Now, Assistant Professor, Division of Engineering Sciences, Purdue 
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Fic. 1. A comparison between the exact solution and several 
approximations obtained by Biot’s method for the melting rate 
in a semi-infinite body under constant heat input. 


To find a minimum of W, we form its first variation and set it 
equal to zero obtaining 


~ 


av = | {H + evT}-sH do (3) 


2 


That a minimum of the function W is found when its variation is 
taken in this way, holding the temperature gradient constant, 
may be seen by forming the second variation 


= f {H+ &VT}-#Hd2+ | (4) 
2 

The first integral will vanish as before, leaving 62% equal to a 

positive definite quantity; hence, the stationary point found 

corresponds to a minimum with respect to the known tempera- 

ture gradient. 


Brot’s METHOD AND A LEAST-SQUARES PROCEDURE 

In Biot’s method, the following procedure is used. A form 
which satisfies all boundary conditions is assumed for the tem- 
perature, say 7), expressed in terms of arbitrary parameters, 
The heat-tlux vector H is then found from Eq. (2) and the arbi- 
trary parameters determined from the differential equations re- 
sulting from satisfying Eq. (3). It is seen that obtaining the 
parameters in this way yields a minimum of W corresponding to 
the assumed temperature, 7), since the temperature gradient, V 7), 
was held constant during the variation. If we assume the exact 
form for the temperature, T,, the above procedure would cause 
W to be a minimum with respect to the exact distribution. How- 
ever, as the temperature gradient which is held constant during 
the variation of V is not in general the same for the two distribu- 
tions, we have no way of comparing the minimums found—i.e., 
there is no means of knowing a priori if the minimum of W cor- 
responding to the exact distribution is any lower than the mini- 
mum of VW corresponding to the approximate distribution. There- 
fore, it is only possible to compare, through the use of the func- 
tion V, two approximate solutioas for the temperature if both dis- 
tributions vield the exact gradient at the instant they are being 
compared. 

Problems arise, however, in which it is necessary to determine 
the best of several approximate solutions at a time when the tem- 
perature gradient is not exactly known. To this end, we intro- 
duce the function ¥ defined by 


v= ae) f {H + AVT}-{H + (5) 
Q 


which is clearly non-negative as it is a perfect square and is 
related to WV by 
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f {aAVT}- {Rv T} do (6) 
Q 


The function Y now provides a means of comparing several 
approximate solutions. When the exact solution occurs—i.e., 


H+ =0 


the function W is identically zero, an absolute minimum. The 
amount by which approximate solutions cause W to differ from 
zero may then serve as a measure of their accuracy, or, as a 
general rule, the smaller , the better the approximation. 

Consider now using the function W in Biot’s method in place 
of yw. Since in forming the variation the method requires that 
the temperature gradient be held constant,* we see that the var- 
iation of ¥ is equal to that of V; thus, the solution obtained will 
be independent of whether we started with either W or Vv. By 
starting with W, however, we see that it is possible to interpret the 
method as requiring a least-squares procedure on the relation H 
+ kVT where the temperature gradient is held constant during 
the variation and energy is conserved in the system. 

As an example, the problem of the melting of a semi-infinite 
body under constant heat input was considered. At the start 
of melting, the function © by itself determines that the best of 
several approximate solutions for a period of time after melting 
has begun will be the one containing the largest number of param- 
eters—i.e. the highest-order approximation. This occurs be- 
cause the set of differential equations resulting from satisfying 
Eq. (3) determine the rate of change of the temperature from 
the initially given or assumed distribution. For this reason the 
solution only holds in some neighborhood of that given or assumed 
distribution. In the example, the temperature is known ex- 
actly at the start of melting, and, therefore, the function 
determines the solution in the neighborhood of this known state 
However, at a later time after melting has progressed, the tem- 
perature will no longer be exactly given by the approximate solu- 
tions; hence, it is not possible to say in general that higher-order 
approximations are necessarily the best. This point was borne 
out by the example which showed that, by the time steady-state 
melting conditions had been achieved, the first approximation 
for the melting rate was better than either the second or third 
approximations. 

Fig. 1 shows the different approximations to the melting rate 
and gives the values of W and wW for all approximations at the 
start of melting and at a later time when steady-state melting 
has been achieved. It is seen from these values that YW cannot be 
used to determine the best melting rate as tr —~ ~, while the use of 
Was r— © clearly shows the first approximation to be the best. 
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* The use of additive functions whose variations are null to form a per- 
fect square, as was done here, arises also in Gauss’s principle of least con- 
straint. See reference 6. 

+ Acomparison with the exact solution‘: for the steady-state melting rate 
was made to determine which was the best solution. 
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First-Order Solution to the Compressible 
Laminar Boundary Layer in Slip Flow? 


Anthony Casaccio 
Senior Research and Development Engineer, 
Republic Aviation Corporation, Farmingdale, L.I., N.Y. 


November 19, 1959 


INCE the phenomena of slip and temperature jump are effects 
S associated with the collisions between gas molecules and the 
body surface, their influence is only important in a region adjacent 
to the body surface, of thickness of the order of one molecular 
mean-free-path length. Therefore, to a first-order approxima- 
tion, the boundary-layer equations of continuum-flow theory are 
still applicable provided the boundary conditions are modified 
to account for the viscous slip and lack of complete accommoda- 
tion at the body surface. 

If the effect of transverse curvature on the boundary layer is 
neglected, the equations governing the laminar motion of a steady 
compressible perfect gas are* 


continuity: (Opur,//Ox) + = O (1) 


momentum: pu(Ou/Ox) + pr(du/dv) = 
— (dp/dx) (2) 


energy: pu(OH/Ox) + pr(OH/dy) = — 
(0/dy) — (3) 


where j = 0, 1 for two-dimensional and axisymmetric flows, 
respectively. 

The modified boundary conditions accounting for first-order 
slip over a constant-temperature surface are given by 


u(x, 0) = T(x, 0) = + (4) 


The coefficients d, and d, depend on the relaxation times neces- 
sary to equilibrate the transfer of momentum and energy and 
on the nature of the reflections which occur between incident gas 
molecules and the body surface. It will be assumed here, to the 
order of approximation considered, that the normal distances 
across which the discontinuities in velocity and temperature occur 
are the same—i.e., d,A = d,A. This latter assumption is not a 
stringent one, since the thickness of the rarefied gas layer ad- 
jacent to the body surface is of order \ so that the coefficients are 
both of order one. Consequently, the boundary conditions to be 
satisfied by the governing equations are 


u(x, 0) = d,(Auy)y-0, v(x, 0) = 0, 
H(x,0) = hy + du(AHy)y=0 
lim u(x, y) = 46, lim H(x, y) = Hs = const. \ 
o 
Introducing a stream function defined by 
Ov /Ov = pury’/po, —(O~/Ox) = pur,’ /po (6) 


and a transformation which combines the Stewartson- and 
Mangler-type variables given by 


x 
0 0 


the governing equations and boundary conditions becomet 


— = SUSU + / 
— = (00/P,) X - (8) 


t+ The author wishes to thank Mr. Robert Sanator for his helpful sugges- 
tions in the writing of this paper. 

* Here rp is the radial cylindrical body coordinate nondimensionalized 
with respect to some characteristic length. 


t Coordinates used as subscripts indicate partial differentiation with ree 


spect to the particular coordinate 
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with 


= O) = pdr/ po) =0 
0) = H(é, 0)/He = 
go + po — O) = 0) = ‘| (9) 


lim yn = U5, lim g = 1 

In these equations, it has been assumed that the Prandtl number, 
P,, is a constant and that the viscosity-temperature variation is 
governed by 

(u/pmo) = C(T/To) 

Since the solution of the governing equations can only be valid toa 
first order in A, perturbation solutions of the form 


v= g= (10) 
are assumed, where ¢ is taken to be (pA,/po). Substituting Eqs. 
(10) into Eqs. (8) and (9) and equating coefficients of «° and « to 
zero, the following two systems of equations result: 


(11) 
with 
| 
Y= = 0, = = 0| 
and 
(12) 


with 
= 

ve = 0, 'g = atn = 0 | 

‘Us, "g = 0 

The zero-order system of equations are exactly those first ob- 
tained by Stewartson and later solved by Cohen and Re- 
shotko.! Moreover, the solution to the first-order system of 
equations may be directly verified to be 


= = gy (13) 


asn— 


Thus, the solutions to the governing equations and boundary 
conditions, taken to first order, are given by 


= + (14) 
g = °g + dur(as/ao)e°g, 

Consistent with these solutions, the local shear and heat flux 
at the body surface can be obtained so that the corrections due 
to slip and temperature jump may be noted. If the heat flow 
due to sliding friction is included in the heat flux from the sur- 
face, there results, to first order, 


— 
(15) 
—q = [(u/Pr)h, + utly=0 = 
(Cyors’Hs/P,)( ps/Poas/ao) =0 (16) 
In these expressions, it should be noted that the correction to the 
local heat flux, due to slip and lack of complete accommodation 
at the surface, is exactly compensated for by the heat flow gener- 
ated by the sliding friction at the body surface. Furthermore, 
if Us is a constant (corresponding to zero pressure gradient ), the 
correction due to first-order slip on the surface shear vanishes, 
yielding the usual flat-plate result. These latter results are in 
agreement with the more recent work of Maslen,? who treats 
the flat-plate slip flow in the presence of an externally induced 
pressure gradient. 
Introducing a similarity parameter, x, such that 


x = W(B/m) f(x) 
U; = °Us~ 


the governing zero-order system of equations become! 


17) 
°g = g(x) 


ite 
lu- 
der 
re 
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be 
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3, 
= 
te a 


fi’ — = (0) (18) 


+ = (1 — P,) X 
My — 1) + Cy — + 9) 


subject to the conditions 


S(O) = f(0) = 0, g(O) = 

lim f’=1, lim g=1f (20) 

& 
The pressure-gradient parameter, 8, is defined as 2m/(m + 1) 
and the velocity ratio, °u/us = °U/Us; = f’, where primes denote 
differentiation with respect to x. Since the right-hand side of 
Eq. (19) is required to be a constant or a function of x—i.e., re- 
quired for similarity, the only cases of practical interest to this 
analysis correspond to P,; = 1 and Mach number factor equal 
to two (locally hypersonic flow). 

For the case P; = 1, the local skin-friction and heat-transfer 

parameters are found to be 


WV Rin: Cro = = 
In In x)[f"(0) — 
V (xm/28)(d In In x)d MoV Ro, (21) 
(Nive 2) = To — = 
V In x)[g'(0)/(ga — (22) 
where Jf, = us/a, and (_), refer to adiabatic wall values. It is 
noted that if the slip-flow parameter 


(MV Rs, 


is zero, these results reduce to the familiar continuum values of 


reference 1. 
For the locally hypersonic case, a more suitable similarity is 
obtained if, (p/p,) ~ ¥," so that 
B = n(1 — y)/y(" + 1) 


Under these conditions, the pressure induced by the displace- 
ment of the boundary layer must be considered. This may be 
accounted for by utilizing an inverse interaction procedure where 
the pressure distribution, p/p, is specified along the edge of the 
boundary layer and the corresponding equivalent body shape 
determined using either the tangent-wedge or tangent-cone 
approximations. The resulting equations governing the local 
skin friction, displacement thickness, and heat transfer, under 
hypersonic approximations are thus given by 


WR wo, Cy, o = */ pat o = 
2niV [(n + 1)/2)(po/p — x(n + 1)/4 X 
(MW (23) 


= 
— 1)/2n4)M ot V + X 


(g — f'2)dx — X 
0 


(N GR 12) = Pr VRot/G = 
ri V + @ — (25) 


x 
0 


and S;, « is the Stanton number. 

If the slip-flow parameter is set equal to zero, Eqs. (23)-(25) 
reduce identically to the corresponding continuum results of 
reference 3 for P, = 1. Owing to the similar solutions approxi- 
mation, the latter considerations are restricted to power-law 
variations of pressure with surface distance. Deviations from 
this type of power law may be treated by piecing solutions 
according to a local-similarity concept. 


where 
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An Additional Note on Equilibrium Skin 
Temperatures 


B. Saelman 

Design Weight Engineer, Lockheed Aircraft Corporation, 
Burbank, Calif. 

December 1, 1959 


| REFERENCE 1, an approximation for the equilibrium skin 

temperature due to a balance between aerodynamic heating 
and radiation was given. It was based on the assumption that 
the terms of higher order than the firstin 7, — 7, [Eq. (21) of ref- 
erence 1] could be neglected. An improved expression results 
when two terms of the binomial expansion are taken. 

Incorporating the term for solar irradiation and assuming that, 
for the body, emissivity equals absorptivity, we have: 


(T, — + = otT.4 = oe[T, — (T, — T.)]4 = 
oe[7,4 — 47,7, — Te) + 67,7, — Te)? .... )~ 
oe(67,77.? — 87,°T. + (1) 
Rearranging terms in powers of 7., 
60eT,?T.? + — 80€T;*)T. + — bnT, = 0 (2) 


Solving for 7, and dividing numerator and denominator by 
oe, we have: 


where 
o= hy /oe 


When only the term involving the first power in 7, — 7, is taken, 
the comparable expression for the equilibrium skin temperature is 


T. = [@7, + (G/o) + + 6) (4) 


As an example, determine the equilibrium skin temperature 
for M = 3.0 at an altitude of 100,000 feet and for a surface mate- 
rial of e = 0.6. For this case, 7, = 390°R., 7; = 1,010°R., and 
hy is approximately 4.02. Neglecting solar irradiation, we have 
by Eq. (4) 


@ = 3.92 X 10° 
Te = 878°R. 
and by Eq. (3) 
T, = 865 °R. 
By trial and error, the exact solution of Eq. (1) is 
T. = 868°R. 


REFERENCE 


1 Saelman, B., Jntegration of Some Thermal Differential Equations, Journal 
of the Aero/Space Sciences, Vol. 26, No. 11, p. 754, November, 1959. 


320 
| 
| 
— | 
; 
<i 
4 
7 
3 
j 
h 
: 
‘ 


ssible 
ACA 


itical 
ining 
sonic 


5-24- 


wdary 
-100, 


REQUIREMENTS for CONTRIBUTIONS 


to the publications of the 


INSTITUTE of the AERONAUTICAL SCIENCES 


The Institute of the Aeronautical Sciences invites both members and nonmembers from any 
country to submit papers for publication in the JouRNAL or THE Agro/Space Sciences and 
Agrro/SPACE ENGINEERING. The Institute, following the practice of other societies, does not pay for 


contributions. 


The following directions for the preparation of papers, if followed by authors, will save corre- 
spondence, avoid the return of papers for changes, minimize the work of preparation for the printer, 
and save the expense due to the charges made for ‘‘author’s corrections,” 


Manuscrirts: The original typewritten copy of the paper is 
desired. Toexpedite review, an additional copy of both the manu- 
script and figures should be submitted. The manuscript must 
be double or triple spaced on one side of white paper sheets, 
consecutively numbered. There should be wide margins to allow 
for the marking of directions to the printer. Correcting, chang- 
ing, or adding to papers after they are in type is costly. It is, 
therefore, imperative that papers submitted be in final form. 
Typographical errors may be corrected on proofs, but if authors 
wish to add material, they may do so at their own expense. In 
mailing, drawings may be rolled, but manuscripts should be sent 
flat. Send by first-class mail (register if you wish for your own 
protection) to the Editorial Office, Institute of the Aeronautical 
Sciences, 2 East 64th Street, New York 21, N.Y. All manuscripts 
will be examined by the Editorial Committee and by the Editor. 
Authors will be advised as promptly as possible whether the paper 
is acceptable for publication. 

Trtizs: The title of the paper should be brief. The name 
and initials of the author should be written as he prefers. The 
use of the full name of an author is advocated because of the fre- 
quent duplication of initials and surnames which sometimes makes 
it difficult to establish the identity of the author. This is par- 
ticularly important for large annual indexing and abstracting 
services, All titles and degrees or honorsare omitted. The name 
of the organization with which the author is associated should 
be placed after his name on a separate line. The date on which 
the paper is received will be inserted by the Editor. The author’s 
title or position should be indicated in a footnote. 

Summarigs oR Apstracts: An abstract to be printed at the 
beginning should accompany each article. It should be adequate 
as an $ and asasummary. It should contain a statement of 
major conclusions reached, since summaries in many cases con- 
stitute the only source of information used in compiling scien- 
tific reference indexes. Abstracts printed in other journals, espe- 
cially foreign, in most cases consist of summaries from printed 
papers. The summary should explain as adequately as possible 
the major conclusions to a nonspecialist in the subject and should 
contain from 100 to 300 words, depending on the length of the 
paper. 

SuBHEADINGS: Subheadings should be inserted by the author 
at frequent intervals. The work of editorial preparation will be 
simplified by the author’s provision of many subheadings. 

Matrer Usuatty DeLetep: Photographs or illustrations of 
little technical interest and not showing advances in general 
practice. Too detailed tabular matter (general results of such 
tables may be included in the text). Lengthy descriptions of 
materials or processes or of preliminary experiments or theories 
that preceded final results; salient features only are of interest. 

REFERENCES: References should be numbered consecutively 
and grouped together at the end of the manuscript, with only 
the corresponding number being mentioned in the text. The 
arrangement should be as follows. For books: ' Durand, W. F., 
Aerodynamic Theory, 1st Ed., Vol. 1, p. 23; Julius Springer, 
Berlin, 1934. For magazines: ‘England, C. R., Crawford, 
A. B., and Mumford, W. W., Some Results of a Study of Ulira- 
Short-Wave Transmission Phenomenon, Proc. IRE, Vol. 20, 
No. 12, pp. 481-482, March, 1933. Please give author, title, 
edition, volume, number, page, publisher, and place and date of 
publication as indicated. Omission of one required fact causes 
much extra editorial work and possible inaccuracies. 


ILLUSTRATIONS: Itlustrations should accompany manuscripts. 
and each should always be referred to in the text by number, 
Drawings or graphs should not be larger than 12 X 16 inches, 
and must be made with jet black India ink on white paper or 
tracing cloth, the latter being preferred. Do not use typewriter 
for lettering. The smallest lettering on 8 X 10 inch figures should 
be no less than '/, inch high. Cross-section paper (white with 
black lines) may be used, but it should not have more than 4 lines 
per inch. If finer ruled paper is used, the major division lines 
should be drawn in with black ink, omitting the finer divisions. 
In the case of finely ruled paper, only blue-lined paper can be 
accepted. Tracing paper and blueprints are not acceptable. 
Lettering and all markings must be large e h to be readable 
after reduction to single-column width (3*/,in.). Mail rolled or 
flat; never fold. Drawings that cannot be reproduced (including 
pencil drawings) will be returned to the author for redrawing, 
thus delaying publication of the paper. Photographs should be 
distinct and show clear black and white contrasts. They must 
be on glossy white paper. Avoid round or oval photographs. 


CAPTIONS AND LEGENDS: Legends or captions must accompany 
each drawing or photograph submitted. If written on the draw- 
ing or photograph, they should be placed below and well outside 
the part to reproduced. Each table should’ have a caption 
such as Table 1, Table 2, Table 3, etc. Captions should be com- 
plete in themselves so as to make the data intelligible to the 
reader without reference to the text. A duplicate list of captions 
for figures should be included as the last page of the manuscript. 
Use ‘‘Fig. 1” (not Figure 1), “Figs. 3 and 4,”’ etc., in both the text 
and the numbering of illustrations. In the text, ““Eq. (1)” or 
“Eqs. (1) and (2)” is used, not “Equation (1).” In captions, 
legends, and in table headings, write all words in full; do not 
abbreviate, except for “‘Fig.”’ and ‘‘Eq.” 


MATHEMATICAL WorE: Formulas may be typewritten or 
carefully written in pen and ink, the writing to be large enough so 
that it can be marked for the printer. Considerable space for mark- 
ing should be allowed above and below all equations. All compli- 
cated equations should be repeated on separate sheets with plenty 
of space left for marking. The solidus should be used for sim- 
ple fractions appearing within the text. Make all expressions 
clear to the typesetter. Greek letters used in formulas should 
be clearly designated by name in the margin of the manu- 
script. The difference between capital and lower-case letters 
should be clearly distinguished and care taken to avoid confu- 
sion between zero (0) and the letter (0), between the numeral 
(one) and the letter (ell) and the prime (’), between alpha and 
a, kappa and k, u and mu, v and nu, n and eta. All sub- 
scripts and exponents should be clearly distinguished. Avoid 
complicated exponents and subscripts. Dots and bars over 
letters or mathematical expressions should also be avoided. 
When it is necessary to repeat a complicated expression, it should 
be represented by some convenient symbol. 


SYMBOLS AND ABBREVIATIONS: The symbols recommended 
in the American Standard Association “Letter Symbols 
for Aeronautical Sciences,” ASA Y10.7—1954, should be used 
wherever practicable. Al! symbols should be clearly written and 
carefully checked. Standard abbreviations should be used, and 
it should be noted that most abbreviations are lower case, such 
as mp.h., b.m.ep., ih.p. b.hp., hp., ete, 
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CoRPORATE MEMBERS 
OF THE 


INSTITUTE OF THE AERONAUTICAL SCIENCES 


AC SPARK PLUG DIVISION, GENERAL MOTORS CORPO- 
RATION 


ACADEMY OF AERONAUTICS 

AEROJET-GENERAL CORPORATION 

AEROLAB DEVELOPMENT COMPANY 

AERONCA MANUFACTURING CORPORATION 

AERONUTRONIC SYSTEMS, INC. 

AEROQUIP CORPORATION 

AGAWAM AIRCRAFT PRODUCTS, INC 

ALLIED RESEARCH ASSOCIATES, INC. 

ALLISON DIVISION, GENERAL MOTORS CORPORATION 

ALUMINUM COMPANY OF AMERICA 

AMERICAN AIRLINES, INC. 

AMERICAN BOSCH ARMA CORPORATION 

AMERICAN STEEL & WIRE DIVISION, UNITED STATES 
STEEL CORPORATION 

AMOCO CHEMICALS CORPORATION 

AMPHENOL-BORG ELECTRONICS CORPORATION 

ANDERSON, GREENWOOD & CO, 

ARDE ASSOCIATES 

ASSOCIATED AVIATION UNDERWRITERS 

AUTOMATION INDUSTRIES, INCORPORATED 

AVCO RESEARCH LABORATORY 

AVIEN, INC, 

BEECH AIRCRAFT CORPORATION 

BELL AIRCRAFT CORPORATION 

BENDIX AVIATION CORPORATION 

BOEING AIRPLANE COMPANY 

BOOZ, ALLEN & HAMILTON 


BULOVA RESEARCH AND DEVELOPMENT LABORA- 
TORIES, INC, 


CANADAIR, LTD, 

THE CESSNA AIRCRAFT COMPANY 

CHANCE VOUGHT AIRCRAFT, INCORPORATED 
THE CHASE MANHATTAN BANK 

CHICAGO AERIAL INDUSTRIES, INC. 

THE CLEVELAND PNEUMATIC INDUSTRIES, INC. 
CONTINENTAL MOTORS CORPORATION 
CORNELL AERONAUTICAL LABORATORY, INC. 
CURTISS-WRIGHT CORPORATION 

DANIEL, MANN, JOHNSON, & MENDENHALL 
THE DECKER CORPORATION 

DEL MAR ENGINEERING LABORATORIES 


DIVISION, GENERAL MOTORS CORPORA- 


DOAK AIRCRAFT COMPANY, INC, 

DOUGLAS AIRCRAFT COMPANY, INC. 

DZUS FASTENER COMPANY, INC, 

EASTERN AIR LINES, INC, 

EATON MANUFACTURING COMPANY 

EDO CORPORATION 

ELASTIC STOP NUT CORPORATION OF AMERICA 
ENGINEERING SUPERVISION COMPANY 
a CAMERA AND INSTRUMENT CORPORA- 


FAIRCHILD ENGINE AND AIRPLANE CORPORATION 
THE GARRETT CORPORATION 

GENERAL APPLIED SCIENCE LABORATORIES, INC. 
GENERAL DYNAMICS CORPORATION 

GENERAL ELECTRIC COMPANY 

GENERAL PRECISION EQUIPMENT CORPORATION 

G. M, GIANNINI & CO,, INC, 

GIANNINI PLASMADYNE CORPORATION 

THE B, F, GOODRICH COMPANY 

GOODYEAR AIRCRAFT CORPORATION 

GRUMMAN AIRCRAFT ENGINEERING CORPORATION 


HYDRO-AIRE, CO, 
INSURANCE COMPANY OF NORTH AMERICA COM. 
PANIES 


INTERNATIONAL BUSINESS MACHINES CORPORATION 

THE INTERNATIONAL NICKEL COMPANY, INC, 

ITT LABORATORIES, DIVISION OF INTERNATIONAL 
TELEPHONE AND TELEGRAPH CORPORATION 


JANITROL AIRCRAFT, A DIVISION OF MIDLAND- 
ROSS CORPORATION 


JOHNS-MANVILLE SALES CORPORATION 

WALTER KIDDE & COMPANY, INC. 

KOLLSMAN INSTRUMENT CORPORATION 

LAVELLE AIRCRAFT CORPORATION 

LEAR, INCORPORATED 

Cc. W. LEMMERMAN, INC, 

LIBRASCOPE DIVISION OF GENERAL PRECISION, INC; 

THE LIQUIDOMETER CORPORATION 

LOCKHEED AIRCRAFT CORPORATION 

LOEWY-HYDROPRESS DIVISION OF BALDWIN-LIMAs 
HAMILTON CORPORATION 

THE MARQUARDT CORPORATION 

THE MARTIN COMPANY 

McDONNELL AIRCRAFT CORPORATION 

MELETRON CORPORATION 

MINNEAPOLIS-HONEYWELL REGULATOR COMPANY 

NATIONAL CREDIT OFFICE, INC. 

NORTH AMERICAN AVIATION, INC, 

NORTHROP CORPORATION 


NUCLEAR DEVELOPMENT CORPORATION OF AMERICA 
PAN AMERICAN WORLD AIRWAYS, INC, 

THE RALPH M. PARSONS COMPANY 

DIVISION, BORG-WARNER CORPORA- 


PHILLIPS PETROLEUM COMPANY 

PLASECKI AIRCRAFT CORPORATION 

RADIO CORPORATION OF AMERICA 
ASTRO-ELECTRONIC PRODUCTS DIVISION 
DEFENSE ELECTRONIC PRODUCTS 


RAMO-WOOLDRIDGE DIVISION, THOMPSON RAMO 
WOOLDRIDGE INC. 


REPUBLIC AVIATION CORPORATION 

ROHR AIRCRAFT CORPORATION 

PAUL ROSENBERG ASSOCIATES 

RYAN AERONAUTICAL COMPANY 

SANDBERG-SERRELL CORPORATION 

SHAFER BEARING DIVISION, CHAIN BELT COMPANY 

SHELL OIL COMPANY 

SIMMONDS AEROCESSORIES, INC, 

SOCONY MOBIL OIL COMPANY, INC. 

SOLAR AIRCRAFT COMPANY 

SOUTHWEST PRODUCTS CO, 

SPACE TECHNOLOGY LABORATORIES, INC. 

R. DIXON SPEAS ASSOCIATES 

SPERRY GYROSCOPE COMPANY DIVISION OF SPERRY 
RAND CORPORATION 

STANDARD OIL COMPANY OF CALIFORNIA 

STANDARD OIL COMPANY (INDIANA) 


STANDARD-THOMSON CORPORATION 
CLIFFORD MANUFACTURING DIVISION 


STANLEY AVIATION CORPORATION 
THIOKOL CHEMICAL CORPORATION 
THOMPSON RAMO WOOLDRIDGE INC, 
TRANS WORLD AIRLINES, INC. 

TURBO PRODUCTS, INC. 

UNION CARBIDE CORPORATION 

UNITED AIR LINES, INC. 

UNITED AIRCRAFT CORPORATION 
UNITED STATES AVIATION UNDERWRITERS, INC, 
VERTOL AIRCRAFT CORPORATION 

VITRO CORPORATION OF AMERICA 
WESTERN GEAR CORPORATION 
WESTINGHOUSE ELECTRIC CORPORATION 
WYMAN-GORDON COMPANY 

YOUNG RADIATOR COMPANY 
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